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1 Introduction 



The issue of verifying the confinement mechanism in Quantum Chromodynamics has 
been a great challenge ever since the phenomenon of string formation between a static 
QQ pair has been conceived by t'Hooft and Mandelstam [IJ] to be a dual Meissner 
effect in the scenario of a Type II superconducting vacuum. Early lattice gauge theory 
attempts to compute the colour field distribution, without recourse to modelling, were 
necessarily limited by the available compute power and lattice methods of the period. 
They rendered qualitative rather than quantitative results, with lattice resolutions, a, 
and quark-antiquark separations, r, restrained to a > .15 fm and r = aR < 1 fm, 
respectively §, |, |, §. 

In recent high precision studies of 577(2) and SU(3) gauge theories || [7| the static 
quark-antiquark potential has been found to be consistent with a linearly rising part 
and a subleading — ir/ (12R) correction as predicted by the bosonic string picture jSL ^ 
for separations above r t ~ .5 fm. Moreover, there are numerical indications for hybrid 



potentials, with gluonic excitations separated by energy gaps nir/R [[H| [TT|] as expected 
from effective string theories ||. 

A compelling evidence about the nature of the confining string from lattice gauge 
theory (LGT) is still lacking. It would require measurements of field distributions 
at quark separations well beyond r t . To study the geometry of the colour flux tube 
between Q and Q sources, one needs an increase both in resolution of the underlying 
lattice and in the linear extent, r, of the string. 

These requirements are not so easily met, since (a) the energy density carries dimension 
a -4 and therefore imposes a lower limit onto the lattice spacing due to statistical noise 
and (b) one is forced to work with very large lattices to attain large quark-antiquark 
separations r = Ra. On top of this, one is of course faced with the ubiquitous problem 
of filtering ground state signals out of an excited state background. 

Thus, in order to really determine the structure of strings in the heavy quark-antiquark 
interaction, one cannot avoid a systematic high precision study, ensuring (a) good 
scaling behaviour as well as (b) sufficient control on finite size effects (FSE) and (c) 
reliable signals for the ground state. 

The superconducting picture for QCD has been modelled in terms of a dual effec- 
tive Langrangian some time ago by Baker and collaborators and worked out subse- 
quently |12[]. Lattice gauge theory in principle offers the laboratory to test such con- 



finement models, as it allows for ab initio studies from the QCD Lagrangian. Within 
the lattice community, there has recently been revived interest to study the role of 
monopole condensation in the confinement mechanism, by recourse to the maximal 
Abelian gauge projection, in SU(2) gauge theory [|13|, [L4], [IB]] . Some first encouraging 
evidence for the dual Meissner effect has also been reported [|T6 1 . Nevertheless, all 



this pioneering lattice work on the confinement mechanism has been carried out either 
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at rather smallish quark-antiquark separations, where the flux tube is not yet really 
developed, or at rather large lattice spacing. 

In fact, there definitely remains a gap: so far flux tubes of sufficient physical lengths 
have never been observed on the lattice! In this paper, we intend to bridge this gap: 
exploiting state-of-the-art lattice techniques — for noise reduction and ground state 
enhancement — as well as the compute power (and memory!) of "small" Connec- 
tion Machines CM-2 and CM-5, we will be able to demonstrate unambiguously that 
quenched SU (2) gauge theory does imply flux tube formation over distances well above 
the 7r Compton wavelength. 

Reliable lattice calculations can only be based on trustworthy error estimates. For this 
reason, we will expose the underlying lattice techniques in quite some detail (Section 
2). There is a shorter Section 3, augmented by three appendices, on (a) weak coupling 
and (b) string model issues that are helpful to appreciate certain qualitative features 
of the field distributions and (c) sum rules for energy and action densities that provide 
an important cross check of the lattice results. The numerical results are presented in 
Section 4 which includes very precise potential data, determination of the Symanzik 
(3 function and many pictures of the flux tubes^. Detailed checks on finite size effects, 
discretization effects and ground state dominance are provided, substantiating the 
interpretation of lattice correlators in terms of continuum fields. Section 5 contains a 
discussion of the shape of the flux tube and the status of Michael's sum rules. 



2 Lattice techniques 

The numerical calculations are performed on lattices with hypercubic geometry and 
periodic boundary conditions in all four directions with volumes L| x ranging from 
16 4 up to 48 3 x 64. Throughout the simulation the standard Wilson action 

Sw = -P £ U^{n) (1) 

n,fi>u 

with 

CV(n) = ^Tr (UMU v {n + £)tfj(n + u)Ut(n)) (2) 
and f3 = 4/g 2 has been used. 

For the updating of the gauge fields a hybrid of heatbath and overrelaxation algo- 
rithms has been implemented |i~7| . The Fabricius-Haan heatbath sweeps [18| have 



been randomly mixed with the overrelaxation step with probability ranging from 1/8 
at (3 = 2.5 up to 1/14 at /3 — 2.74. The links have been visited in lexicographical 
ordering within 2 4 hypercubes, i.e. within each such hypercube, first all links pointing 
into direction 1 are visited site by site, then all links in direction 2 etc.. 

3 We regret being unable to expose the colour flux tube in colour in this medium! A data base of 
colour images is under construction. It can be accessed via anonymous ftp from wptsO.physik.uni- 
wuppertal.de. The (compressed) .rgb and .ps files are deposited in the directory pub/colorflux. 
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2.1 Prerequisites 



In order to substantiate continuum results from lattice calculations it is of utmost 
importance to investigate the impact of the finite lattice volume as well as the scal- 
ing behaviour with the lattice spacing a. This requires simulations both at (a) fixed 
lattice coupling (3 (i.e. spacing a) with a varying number of lattice sites and (b) (ap- 
proximately) fixed physical volume but different lattice resolutions. 

Of course one wants to work — within the computational means — as close as possible 
to the continuum limit, i.e. at as large /3-values as possible. The bottleneck is set by 
the memory requirements due to the increase of the number of lattice sites (needed to 
compensate for a smaller a) as well as by the computer time required to suppress the 
statistical noise. Since the operators under investigation scale with the fourth power 
of the lattice resolution (up to In a terms from anomalous dimensions, see below), the 
latter limitation is the more serious one, restricting all preceding lattice studies to 
P < 2.5. 

Though simulations at small (3 values allow for rather large physical volumes, lattice 
artefacts are expected to spoil results at physically interesting scales. Moreover, we 
should mention that our smearing procedure provides inferior ground state overlaps 
at large lattice spacings. There is more reason to stay away from too coarse lattices: 
one needs a sufficiently large T-range for verification of T-plateaus in the bona fide 
physical quantities. In addition, at smaller values of j3, the lower limit T > 3, implied 
by the minimal temporal extent of the (□)w operator, amounts to overly large physical 
separations and leads to small signals of the Wilson loops. 





P = 2.50 


P = 2.50 


P = 2.635 


P = 2.74 




16 4 


32 4 


48 3 x 64 


32 4 


K 


.0350(4) 


.0350(12) 


.01458(8) 


.00830(6) 


a/fm 


.0826(5) 


.0826(14) 


.0541(2) 


.0408(2) 


a"7GeV 


2.35 (1) 


2.35(4) 


3.64(1) 


4.83(2) 


aLs /fm 


1.32(1) 


2.64(4) 


2.60(1) 


1.31(1) 


MC sweeps 


868000 


255600 


53800 


152000 


thermalization sweeps 


2000 


51000 


4200 


4000 


meas., sweeps between meas. 










Wilson loops 


8680, 100 


2046, 100 


248, 200 


1480, 100 


colour flux distribution 


8680, 100 


2046, 100 


248, 200 


670, 100 


plaquette 


86800, 10 


20460, 10 


992, 50 


14800, 10 



Table 1: Simulation parameters. The physical scales have been computed from the 
value \f~Ka = 440 MeV. 

In short, one has to find a way between Scylla and Charybdis: i.e. to compromise 
between the shortcomings of both small and large lattice spacings a. We have chosen 
to simulate at P — 2.5,2.635, and 2.74 at various lattice volumes, ranging up to the 
unpreceded volume 48 3 x 64. Our simulation parameters are summarized in Tab. [I]. 
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As a by-product we compute the static potential and obtain the most precise set of 
SU (2) string tension values that has ever been computed on the lattice. Details of the 
fitting procedure are explained in Section [4.1. 2| . The (lattice) string tension K relates 



the lattice spacing a to a physical scale: k = Ka 2 . We ascribe the "canonical" value 
y/K = 440 MeV to the square root of the string tension. Needless to say that this 
scale, taken from real world QCD Regge trajectories, only serves as an orientation for 
its poor man's quenched two colour version: SU{2) gauge theory. Nonetheless, the 
quantitative agreement between the SU(2) and SU(3) potentials is remarkable. We 
also point out that the effective string model with which we are going to compare our 
results does not depend on the underlying gauge group. 

From the string tension measurement we find the following (approximate) relation 
between the present lattice spacings: a^ 5 : a^g 35 ■ a 2.7A ~ 1 : 1-5 : 2. Thus, the 16 4 
lattice at f3 — 2.5 has approximately the same physical volume as the 32 4 lattice at 
(3 = 2.74. The same holds true for the 32 4 lattice at (3 = 2.5 and the 48 3 x 64 lattice at 
(3 = 2.635. These pairs of lattices can be used to investigate the a (in) dependence of 
the results. In order to reveal possible volume effects, the 16 4 and 32 4 lattice outcomes 
at /3 = 2.5 will be compared with each other. 

As a prerequisite to the present investigation, let us consider the static QQ potential 
which can be computed from Wilson loops, W(R,T). A Wilson loop, i.e. an ordered 
product of link variables along a closed rectangular path with spatial separation R and 
temporal extent T, can be interpreted as a world sheet of a QQ pair: at Euclidean 
time r = a creation operator 

T R = Q(0)C/(0 -> R)Q^(R) (3) 

with a gauge covariant transporter U(0 — > R) is applied to the vacuum state |0). 
The QQ pair is then propagated to r = T by static Wilson lines in presence of 
the gauge field background, and finally annihilated by application of Tr. A spectral 
decomposition of the Wilson loop exhibits the following behaviour (T = ] denotes 
the transfer matrix, T\n) = e~ En \n))\ 

Tr (T R T T T R T L ^ 



(W(R,T)) 



Tr (T Lt ) 

L-^E \(m\T R \n,R)\ 2 e- v ^ T e- E -^ (4) 



■dm ^ m,n 



= J2 \d n (R)\ 2 e- y ^ R)T x (l + O (e-^ 1 ^ 

n 

where d n (R) = (0\T R \n, R). \n,R) is the nth eigenstate in the charged sector of the 
Hilbert space with non-vanishing overlap to the creation operator r^, while \n) is the 
nth eigenstate of the zero charge sector. V n (R) denotes the nth excitation of the QQ 
potential and the vacuum energy Eq has been set to zero. E% is the mass gap, i.e. the 
mass of the Af glueball. 

Actually, we are not restricted to on-axis QQ separations, R = (R, 0, 0). Planar 
Wilson loops can be easily generalized to off-axis separations by connecting sources 
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that do not share a common lattice axis. In the present investigation, the following 
off-axis directions have been realized: 

di = (1,0,0) , d 2 = (1,1,0) , dg = (2, 1,0) , 

d 4 = (1,1,1) , d 5 = (2,l,l) , d 6 = (2,2,l) , (5) 

with separations m^dj up to 7711,7712,7714 < Ls/2 and 7723, 7725, m6 < vLs/4. For the 
largest lattice, L$ = 48, this amounts to a measurement over a set of 108 different 
separations. All paths have been chosen as close to the shortest linear connection 
between the sources as the lattice permitted. 



2.2 Noise reduction 



In this section we will shortly discuss the implications of noise reduction that we 
achieved by integrating out the temporal links in the Wilson loops analytically |19| . 

The link integration amounts to the following substitution: 

fcrr^dU Ue? s ">^ 

u,{n) - v ' (n) = (6) 



with 



S n jU) = jTr (UFl(n)) , F,{n) = Y.U u {n)U,{n + 0)Ul{n + (i) . (7) 



2 

V^n) is in general not an SU(2) element anymore 



In this way, time-like links are replaced by the mean field value they take in the neigh- 
bourhood of (stroboscopically frozen) links that interact through the staples F^n). 
Only those links that do not share a common plaquette, can be integrated indepen- 
dently. This holds in particular for all temporal links within our spatially smeared 
Wilson loops, iff R > 1. 

In case of 577(2) gauge theory, V±(n) can be calculated analytically: 
I n denote the modified Bessel functions. 

The statistical error, AO, of an observable, (O), calculated without link integration, 
is related by a constant s < 1 to the corresponding error with link integration, AOu = 
sAO. In order to discuss the impact of link integration on noise reduction, we start 
from the naive expectation that each integrated link contributes equally to s, i.e. 
we assume s = x 2T with 2T being the number of integrated links used within the 
construction of (O). 
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In order to estimate the value of s, let us consider on-axis Wilson loops with inte- 
grated temporal links. On symmetric lattices (Lt = Lg) we expect from the relation 
(W(R,T)) = (W(T,R)): 

^W(R,T) {T _ R) 

AW(T,R) ' K} 



This leads to the estimate for x . 



1 (AW(R,T) 



^ = exp|— .og^^^M . (10) 



Our data (with bootstrapped errors of the errors) is consistent with a factoiQ x = 
.889 ± .001 (x = .890 ± .001) for (3 = 2.5 (J3 = 2.635). Thus, application of link 
integration at time T = 6 (the largest temporal extent, used in the computation of 
the colour field operators, see below) amounts to a reduction of computer time by a 
factor x" 24 = 16.8 ± 0.4. 



The improvement achieved by link integration tends to be smaller at smaller lattice 
spacings. This is due to the fact that the physical extent of the neighbourhood to be 
integrated out becomes smaller. On the other hand, the error of a non-link integrated 
operator, measured on lattices with constant physical volumes but different couplings, 
also decreases with the lattice spacing (temperature (3~ l ). At the bottom line, the 
two effects almost cancel each other and the relative errors of link integrated Wilson 
loops appear to remain rather independent of the lattice resolution, provided that the 
physical lattice volumes and the number of measurements are kept constant. 



2.3 Ground state enhancement 



The physically interesting ground state potential, V(R) = Vq(R), can be retrieved in 
the limit of large T: 

(W(R,T)) = Y / C n (R)e- y ^ T ^C (R)e- y ^ T (T - oo) . (11) 

n 

The overlaps C n (R) = |<i n (-R)| 2 > obey the following normalization condition: 

£tf„(i2) = l ■ (12) 

n 

The path of the transporter [7(0 — > R) used for the construction of the QQ creation 
operator (Eq. ^) does not affect the eigenvalues of the transfer matrix and is by no 
means unique. One can exploit this freedom to maximize the ground state overlap by 
a suitable superposition of such paths, aiming at Cq(R) ~ 1. At any given value of 
R, the final deviation of Co(R) from one can serve as a monitor for the suppression of 
excited state contributions actually achieved in this way. 

4 For small R and T where the statistical errors have reached the same order of magnitude as 
the numerical machine accuracy, we performed additional computations of AO and AOu on smaller 
subsamples to ensure that the errors still follow the statistical 1 / \f N meas expectation. 
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Figure 1: Visualization of a smearing iteration. 

In the present simulation an iterative procedure (with nj ter iteration steps) has been 
applied [ 
a "fat" link, 



2jJ : each spatial link U{(n), occurring in the transporter, is substituted by 



Ui(n) — > N [ a^(n) + £ ^(n)C^(n + + 



(13) 



with the appropriate normalization constant AT and free parameter a. One such it- 
eration step is visualized in Fig. |l[ For this smearing, the links are visited in the 
lexicographical ordering of the updating sweep. 



We find satisfactory ground state enhancement with the parameter choice nu er 
and a = 2. 
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One can define approximants to the asymptotic potential values and overlaps, V(R, T) - 
V(R) and C (R,T) — > C (R) (T — > oo). Due to the positivity of the transfer matrix 
T, these quantities decrease monotonically (in T) to their asymptotic limits: 



(W(R,T)) T+1 



Ci(i?) 



h{R,T) + 



C (R,T) 



(W(R,T + 1)) T 



C (R) 

C (R) +Ci(R)h(R,T) + . 



with 



h(R, T) = e~ Av ^ T (l - 



-AV(R) 



AV(R) = Vl(R) - V(R) 



(14) 
(15) 

(16) 



In our analysis, we follow these approximants until they reach a plateau. 



As we wish to maximize Cq(R), we would like to acquire a qualitative understanding 
of the underlying physics. For this purpose, we consider unsmeared on-axis Wilson 
loops. Combining Eq. with the R-T symmetry^ (W(R,T)) = (W(T,R)) and the 
parametrization of the potential V(R) = Vq — e/R + KR we obtain 



ln(W(R,T)) w ]nC {R,T)-V T + eT/R- KRT 

= In D - V (R + T) + e(R/T + T/R) - KRT 



(17) 



J This symmetry is only exact on lattices with L$ = Lt- However, within statistical accuracy it 
also holds true on the 48 3 x 64 lattice. 
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Figure 2: TTie ground state overlaps Cq(R) versus R at (3 = 2.635. In addition to 
the overlaps of smeared operators (diamonds), finite T approximants, C (R,T) for the 
overlaps from unsmeared (on-axis) operators are plotted. The dashed line denotes the 
extrapolated large T limit for the unsmeared overlaps and should be a valid approxima- 
tion to the large R behaviour. 

for large R and T and arrive at the estimate (with constants Vo and e obtained from 
the potential analysis below), 

C (R,T) = De-( v °- e/T)R or C (R) = De~ v ° R . (18) 

This parametrization of the unsmeared overlaps turns out to describe the off-axis data 
too if we allow for a smaller (direction dependent) constant Vq. 

The self-energy term Vq/o diverges in the continuum limit, a — > 0. Thus, the overlaps 
at fixed physical separation^], r = Ra, decrease with increasing (3. This feature is in 
accord with the following consideration: in the scaling region the transverse size of the 
QQ wave function is expected to remain constant in physical units while the transverse 
extent of the string like creation operator remains on the scale of the lattice resolution. 
Thus, the ground state overlap of this operator decreases with increasing correlation 
length. 

The ground state overlaps of smeared Wilson loops at (3 — 2.635 are shown in Fig. 0, 
together with on-axis approximants to the unsmeared overlaps and the asymptotic 
(large R) estimate of Eq. [18] (dashed curve), where the coefficient D pa 2.3 was 
obtained by fitting the large T data to Eq. |18|. There is obviously a dramatic im- 

6 At fixed (lattice) R, a (slight) increase is observed and expected. 
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provement from the use of an extended creation operator. In case of the smeared 
links, the (unsmeared) self-energy contribution Vq appears to be reduced to a number 
/ C V that is sufficiently small to allow for an expansion of the exponential factor 
Co(R) = D exp(-fR) m D(l—fR): the ground state overlaps of smeared Wilson loops 
exhibit a linear R dependence throughout the observed R region. Moreover, rotational 
invariance in terms of the overlaps is restored for all on- and off-axis separations. 

For P = 2.5 the overlaps vary between C (V2a) = .95 and C (r m ) = .73 on the 32 4 
lattice and from Co(a) = .98 to Co(r m ) = .81 on the 16 4 lattice. Within the same 
physical region the (3 = 2.635 overlaps range from Co(V2a) = .98 to C (r m ) = .81. At 
(3 = 2.74 we have used an inferior set of smearing parameters {n iter = 40 and a — 0); 
yet we achieve overlaps of C (a) = .96 and C(r m ) = .84. We have chosen r m ps 1.2 
fm for the comparison. This scale corresponds to r m a _1 = 8\/3, 12\/3, 16a/3 for the 
three /5-values, respectively. Even at fixed physical r the overlaps tend to increase 
with /3, unlike in the situation with unsmeared operators: the wave function becomes 
smoother at increased correlation length and can be better modelled by the iterative 
smearing procedure. For the largest distance realized (2.25 fm at f3 — 2.635) we still 
achieve the value C (24-\/3a) = .72. 

The success of smearing is twofold: (a) for rather small values of T, extraction of the 
ground state potential becomes possible and (b) the signal-to-noise ratio is greatly 
improved as C (R) (and the signal) increases, especially for large values of R. 



2.4 Lattice determination of colour fields 

The central observables in our present investigation are the action and energy densities 
in presence of two static quark sources (with separation R) in the ground state of the 
binding problem: 

^(n) = ±(£ R (n)+B R {n)) , (19) 
Mn) = \ (£ R (n) - B fl (n)) , (20) 

with 

S R {n) = (E 2 (n))| , RH o> , B R (n) = (B 2 (n))| 0ii?Ho) (21) 

and 

(O)|o,fl>-|o> = (0, R\O\0, R) - (0|O|0) . (22) 

The sign convention corresponds to the Minkowski notation with metric 
r] = diag(l, — 1, — 1, — 1), in which B R (n) < 0, £ R (n) > 0. We point out that the 
Minkowski action density carries a different sign relative to the (negative) Euclidean 
action, i.e. S w = - £ n |(E 2 - B 2 )| 0) . 

We shall extract these observables from the correlations between smeared Wilson loops, 
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W = W(R,T), and (unsmeared) plaquettes D(r) = C/^(n,r) (Eq. |) 7 , 

(oH»4 WD|T/H ;; D(r/2 - sl)l -(o) . (23, 

5" denotes the distance of the plaquette from the central time slice of the Wilson loop 
and takes the values S = 0, 1, . . . (S = 1/2, 3/2, . . .) for even (odd) T. 

The plaquette insertion acts as the chromo dynamical analogue of a Hall detector 
in electro dynamics^. For < S < T/2, (□(5 , ))w can be decomposed into mass 
eigenstates as follows 



Tr [T [rT/2+S u j-T/2-S + j-T/2~S D j-T/2+S^ r f T L T ~Tj 

(□(S)) w = 2Tr (rr r rtT^-T) ( D ) 

= (0,i2|D|0,i2) - (0|D|0> 

+ 2Re ^(I,i2|n|0,i2)^ e~ AVT/2 cosh(AVS) 

+ ^((l,R\n\l,R)-(Q,R\n\0,R))e- AVT (24) 

I dor 

+ 2Re (^-(2,R\n\0,R)^j e "»-W J cosh((y 2 - 7)5) 
+ 0( e - AV -(l T - s )) . 

A 7 denotes the gap between the ground state and the first excitation (Eq. |16|). In 
principle, \d n \ 2 = C n (R) and V n (R) can be determined from smeared Wilson loops. 
The non-diagonal 0(e~^ Vn ~ v ^ 1 V 2-5 )) coefficients can only be obtained from a fit to 
the time dependence of the above operator. As we shall see in the next paragraph, 
(y 2 — V)/2m AV. Thus, a measurement of the excited state colour field distribution 
(□)|i,i?)-|o) appears to be unfeasible with the present method. 

In string model calculations the separations between ground and excited state poten- 
tials without gluonic angular momentum, i.e. within the A\ g representation of the cubic 
symmetry group, are found to be multiples of 2tt/R [|8], |TD[. This feature is in accord 



with numerical simulations of SU(2) and SU(3) gauge theories [I0[ [LI]. Therefore, as 
a net result, we expect the following asymptotic behaviour: 

(n(s)) w = (n) lo , RH o) 



7 We do not follow the authors of Ref. Q] who, in order to reduce statistical fluctuations, advocate 
to subtract (Wd(n)) / (W) with the reference point, n, taken far away from the sources rather than the 
vacuum plaquette expectation (□). In this way, we avoid possible shifts of the normalization relative 
to the vacuum energy and action densities. We would like to point out that we found no reduction 
in statistical errors for smeared Wilson loop operators by using the above suggestion. However, we 
have been able to confirm this observation of Ref. B for unsmeared Wilson loops. 

8 We note in passing, that the authors of Ref. E2] have chosen to connect the plaquette to the 
Wilson loop via two Wilson lines and take one overall trace instead of two separate ones, as this leads 
to an improved signal to noise ratio. However, a prove that this observable indeed can be interpreted 
as a colour field density in presence of a static QQ pair is missing. Moreover, the constraint through 
sum rules is lost. 



10 



+ Cl e~ nT/R cosh{27c S/R) 

+ c 2 e- 2wT/R (1 + c 3 cosh(4vr5/ J R)) + 



(25) 



with q being free parameters. They are, contrary to the coefficients of the spectral 
decomposition of the Wilson loop (Eq. [TT]) , not necessarily positive! Be aware, that 
Ci vary with the QQ separation as well as with the spatial position n of the plaquette 
insertion. 

The deviations from the asymptotic values are governed by 0(e~ Av ^ T ^ 2 ~ s ^) terms, 
compared to order e~~ AVT corrections in case of the potential (Eqs. [14], [161). So, the 
issue of optimization for ground state dominance is certainly more critical for field 
measurements. While the reduction of systematic errors would ask for large T-values, 
the suppression of statistical uncertainties would lead one to the contrary. Obviously, 
the reasonable strategy is to ensure that systematic and statistical errors are kept in 
balance. 



A weak coupling expansion of the plaquette yields the square of the Maxwell field 
strength tensor = F^ v a c /2: 

U, u = I - + 0(a 6 ) . (26) 

Thus, by an appropriate choice of □ = Ua (□ = Uj k ) expectation values of squared 
chromo electric (magnetic) field components can be obtainedP], in the limit of large T: 

%(Uu(S)) w T ^ (^(n))| 0j * H o) , (27) 

- 2 -^(U 3k {S)) w ^ (^ 2 (n)> M Ho) , M = l . (28) 

The finite T corrections to these relations have been elaborated in Eq. |2f| Note, that 
Ef = EfEf = 2Ti Ef. 



2.5 Implementation of colour field operators 

For measurement of the colour field distributions, the appropriate plaquette operators 
are suitably averaged in order to obtain chromo magnetic or electric insertions in 
symmetric position to a given lattice site, n. For the electric insertions two plaquettes 
are averaged: 

C4i (n) -> l - (C/ i4 (n - eO + C/ i4 (n)) . (29) 
For the magnetic fields four adjacent plaquettes are combined: 

U jk {n) -> - {U jk (n - ej - e k ) + U jk (n - e,-) + U jk (n - e k ) + U jk (n)) . (30) 

9 Remember, that we do not obtain any information on the components of E and B themselves 
since, in general, (O 2 ) ^ (O) 2 - 
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Notice that while Bf is measured at integer values of r, Ef is measured between two 
time slices. To minimize contaminations from excited states (Eq. |25|), r is chosen as 
close as possible to T/2. For even temporal extent of the Wilson loop, this means 
S = for the magnetic field operator and S = 1/2 for the electric field insertion, while 
for odd T, S = 0(1/2) for electric (magnetic) fields. 

For measurement of the colour field distributions we have restricted ourselves to on- 
axis separations of the two sources. All even distances R = 2, 4, . . . , -R max with -R max = 
8, 24, 36 for Ls = 16, 32, 48, respectively, have been realized. In order to identify 
the asymptotic plateau, T was varied up to0 T = 6. The colour field distributions 
have been measured up to a transverse distance n± = 2 along the entire QQ axis. 
In between the two sources and up to 2 lattice spacings outside the sources, the 
transverse distance was increased to n±^ max = 6, 10, 15 for the three lattice extents 
Ls = 16,32,48, respectively. In addition to "on-axis" positions, n = (rii,n 2 ,0), we 
chose plane-diagonal points n = (ni,n 2 ,n 2 ) with n 2 < n± >max /\^2. We averaged 
over various coordinates n, exploiting the cylindrical and reflection symmetry of the 
problem. All this amounts to 9576 (26244) different combinations of R, T, and n on 
the 32 4 (48 3 x 64) lattices, on which both (Ujk)w an d (Ui^w have been computed. 

The temporal parts of the Wilson loops, appearing in the colour field correlator, have 
been link integrated. Therefore, the electric components have only been determined 
at distances larger than one lattice spacing away from the sources. For the case of the 
32 4 lattice at /3 — 2.5 we have substituted the missing values by the corresponding 
entries, computed on a 16 4 lattice without link integration. Distances so close to the 
sources are not relevant to continuum physics anyway, due to contamination from the 
heavy quark self-energy and lattice artefacts. 



3 Theoretical expectations 



3.1 Perturbative scenario 



A perturbative order g A computation of the lattice potential can be found in Ref. [23 
Here, we recall the one gluon exchange result only: 

V(R) = -C F g 2 (G L (R) - G L (0)) 2=3 -C F g 2 -±- (31) 

AirR 

where we have dropped the (divergent) self energy in the continuum expression. The 
lattice gluon propagator G L (R) (Eqs. [78], [77]) can be computed numerically on finite 
lattices. For SU{2) one has C F = (N 2 - l)/(2N) = 3/4. For completeness, this 
expression is derived in Appendix 0. A renormalization of the bare lattice coupling 

10 For historical reasons, on our small lattice volumes (16 4 at j3 = 2.5 and 32 4 at j3 = 2.74) only the 
odd values T = 1, 3, 5, 7 have been realized. 
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g 2 = 2N/ (3 turns out to be the main effect of the loop diagrams that occur in the next 
order. 



In order to investigate the nature of lattice artefacts, we have performed a weak cou- 
pling expansion of the action and energy densities. The lowest order term is a two 
gluon exchange. To this order the action and energy densities turn out to identical, 
both receiving just contributions from electric plaquettes. Details of the calculation 
are contained in Appendix The lattice integrals of the result (Eqs. [53], [S3]) have 
been computed numerically. 
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Figure 3: Comparison between continuum and lattice dipole fields in the center plane 
between two sources, separated by R = 12. The ordinate, n, is the distance from the 
QQ axis. The values have been multiplied by a factor (47r) 2 /{Cf9 2 )- Crosses and the 
solid line correspond to the infinite volume results. Squares and the dashed line are 
obtained at finite volume, R/Lg = 12/32. 

In the continuum limit one finds the expression (from Eq. RO), 

e£\0,n ± )=g 2 C F , A l, 2/rnni f i 2 x3 , (32) 



(4tt) 2 (B?/4 + 



n 



for the energy density distribution in the central transverse plane. In Figs. |3] and f| we 
present a comparison of the dipole fields on finite (32 3 ) and infinite lattices with their 
continuum forms^], for separations R = 12 and R = 4, respectively. The field positions 
are chosen both along a transverse lattice axis and a plane-diagonal (multiples of a/2). 



ii 



The continuum formula on a finite lattice has been elaborated in Appendix 
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Figure 4: Same as Fig. [| for R = 4. 



Up to order g 4 corrections, perturbative lattice and continuum calculations equally 



As argued in the appendix, perturbation theory is expected to describe the energy 
density better than the action density. 



3.2 Nonperturbative expectation 

In the limit of large QQ separations, i.e. if the width of the flux tube becomes small 
relative to its length, an effective relativistic string model is expected to describe the 
infra red aspects of the interaction. Classical solutions of such string Lagrangians 
predict, in agreement with the strong coupling expectation of pure gauge theory, an 
area law of Wilson loops and, thus, a linearly rising long range contribution to the 
potential. However, in reality, a quantum mechanical string will fluctuate. An ultra 
violet cut-off has to be imposed on the wave length of such fluctuations, beyond which 
longitudinal degrees of freedom become important and the (non renormalizable) string 
theory looses its applicability. For a huge class of string models the string fluctuations 
lead to a universal subleading Coulomb type contribution ||, — (d — 2)7r/(24i?), to 
the potential in the Gaussian approximation (d denotes the number of space-time 
dimensions). For large R, excitations are expected to be separated from the ground 

12 To obtain the continuum expression, X)n fl3 simply has to be replaced by J d 3 x. 



lead to0 




(33) 
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state by multiples of tt/R |§. 

The leading order expectation of string models for correlators of smooth, large Wilson 
loops, W and □ with boundaries <9W and do is 

(WD) - (W) (□) oc exp (-KA(dW, 3D)) , (34) 

where A(dW, <9d) is the minimal area of a surface with boundary dW U do and .K" is 
the string tension. Approximating the Wilson loop and the plaquette by circles, the 
authors of Ref. obtained in the limit n± -C R 

( E D « ex P (-^l) ■ ( 35 ) 



lni? 

Thus, the central width of the fluctuating string 

Jdn± n 3 ± e R (0,n± 
fdn± n±e R (0,n ± ) 

is expected to diverge logarithmically with the quark separation: 



= ( n ±>^ = FX: I ~ \ ( 36 ) 



C = 5 oln^ , (37) 

where R® is an ultra violet cut-off parameter. In a quantum mechanical calculation, 
this relation has been confirmed in the Gaussian approximation for the probability 
of the fluctuating string, crossing the central transverse plane n\ — at the position 
n± H3- 



For small distances, the perturbative result of Eq. [32] suggests a linear divergence of 
the width: 

For rij_ R (and large r) where the string picture is applicable a Gaussian transverse 
profile of the flux tube is expected. At large n±, however, correlators of (unsmeared) 
Wilson loops with plaquettes can be viewed as glueball correlation functions in rotated 
space-time. Thus, for n± ^> R,T an exponential form, governed by the mass gap, 
might be expected^- 



3.3 Sum rules 



Some important consistency conditions, relating the local chromo field operators to 



the global QQ potential have been derived by Michael [25|. In the following we will 



13 However, the wave function of the QQ pair, created at r — has to be decayed into its ground 
state, before the colour fields are measured at r = T/2. Due to the structure of the action, a 
Hamiltonian evolution in the strong coupling limit only allows hopping between neighbouring sites. 
Thus, in this limit, communication occurs only between sites within the light cone n±_ < T/2 and the 
limit n_L S> T is not justified. As illustrated by the above example, the exponential decay prediction 
for large n± has to be taken with care. 
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shortly recall these sum rules. More details and comments related to contaminations 
from excited states can be found in Appendices |B] and |Cj 



The action sum rule relates the action to the derivative of the potential with respect 
to the inverse coupling 

Eo v H(n) = _i<m (39) 

n a amp 

9lna AR)--^ , (40) 



d In (3 a din (3 

where we have decomposed the potential V(R) = av(R) + Vo into a physical part, v(R), 
that remains constant as a — > 0, and a diverging self energy contribution, Vq. The 
action sum rule is an exact relation. It is in accord with the perturbative expectation 
Eq. |33], which follows by inserting the leading order expression for V(R) (Eq. |31]) into 

The energy sum rule involves derivatives with respect to anisotropic spatial and tem- 
poral couplings. After relating the latter to the isotropic lattice coupling, /3, via a 
weak coupling series, one ends up with an approximate sum rule. Thus, unlike the 
action sum rule, the energy sum rule is not exact on the lattice. Here, we just state 
the leading order expectation 

E = + ^) + 0(P x )) . (41) 

The correction to this energy conserving rule reflects the fact that the local plaquette 
operator undergoes a renormalization. However, mean field arguments (Appendix |A|) 
suggest only small corrections. The energy sum rule is also in accord with Eq. |33[ 

The leading order contribution to the self energy V , CpGL(0)g 2 , merely changes sign 
when differentiated with respect to In (3. As a consequence, this contribution to both 
the action and the energy sums diverges like l/(alna). Due to the localization of the 
self energy to the vicinity of the two sources, the peaks of the distributions diverge 
like l/(a 4 lna) in physical units (or like 1/lna when measured in lattice units). The 
physical part v(R) on the r.h.s. of the action sum rule is accompanied by an anomalous 
dimension d In a/d In (3 oc In a. For this reason, the measured lattice action density aa A 
is expected to scale like a 4 In a outside the peaks while the energy density vanishes like 
a 4 (in lattice units). £ and B mix under renormalization group transformations since 
the sum of both densities carries dimension a 4 while its difference is accompanied by 
a 4 In a. Thus, only the energy density and the combination {d In (3/d In a)a are relevant 
to the continuum limit. 



4 Results 



16 



4.1 Static potential 
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Figure 5: V(R, T) as a function of T for R = 20, R = 10, R = 5 and R = 2, 
respectively, at (3 = 2.635. The plateau values obtained at T = 5 are indicated by the 
horizontal dashed lines. 

To prepare the stage for the colour flux investigations, we had to calculate the static 
QQ potential. It will render the string tension, which serves to fix the physical scale 
and relates the results, obtained at different (3 values to each other. 



4.1.1 Data 



The potential data has been obtained by the method described in Section [2.3[ Finite T 
approximants V(R, T) and Cq(R, T) to the ground state potential values and overlaps 



are computed according to Eqs. [13] and 15. These are traced until a plateau (in T) 



is reached. The numerical situation is illustrated in Fig. |5] for a few typical quark 
separations at (3 = 2.635. For the 16 4 lattice at (3 — 2.5 the T min = 3 approximant 
has been found to agree with the plateau values while for the 32 4 lattices at (3 = 2.5 
and (3 = 2.74, T min = 4 had to be taken and at /3 — 2.635 we went as far as T min = 5. 
To exclude any remaining systematic bias on the fitted parameters, all fits have been 
performed for T = T min , and T = T min + 1. Within statistical errors and fixed R range, 
the fit parameters remained stable. For larger (3 values, the reduced physical t = Ta 
separations appear to be partly compensated by better ground state overlaps. 



We note, that the actual value of T m i n is not only affected by the ground state over- 
laps but also influenced by statistical errors that depend on the particular number of 
measurements, physical volume and method (link integration). 

In Tabs. Q-f| we have collected results on the potential values, V(R), and overlaps, 
Cq(R), up to a physical distance of about .7 fm. This scale has been obtained from the 
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beta=2.50 , Ls=16 h©— i - 
beta=2.50 , Ls=32 h — i 
beta=2.635, Ls=48 hb— i 
beta=2.74, Ls=32 h^-h 



_3 lj i i i i i i i 

0.5 1 1.5 2 2.5 3 

R sqrt(K) 

Figure 6: The potential, measured on the four lattices, scaled in units of the string 
tension. The solid line refers to the string picture expectation V(R) = K R — ir / (12R) . 
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Table 2: Potential and overlap values at (3 = 2.5. 
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Table 3: Potential and overlap values at (3 = 2.635. 

relation k = yKa = 440 MeV. For larger separations we refer to the parametrizations 
presented below since no systematic deviations are observed from the interpolating 
curve (that is dominated by the linear part of the potential). Remember, that all esti- 
mates for the potential and overlaps constitute strict upper limits to their asymptotic 
(T oo) values. The paths, displayed in the second column, are numbered according 
to Eq. |. 

In Fig. |^ we show the familiar scaling plot for the potential in form of the combina- 
tion (y(RyK) — Vo) / y/K with Vq and K as obtained from the four- parameter fits, 
described below. Notice, that we can trace the potential up to the impressively large 
separation of 2.3 fm! The curve represents the string picture prediction R — The 
nice scaling between the potentials illustrates the restoration of continuum rotational 
invariance at remarkably small lattice separations. 

4.1.2 Potential fits 

Our potential values have been fitted to the parametrizations 

V(R) = V + KR-^ (42) 
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Table 4: Potential and overlap values at (3 = 2.74. 
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(44) 



where Gl(R) is the lattice gluon propagator (Eqs. [7^, [77|), computed on an infinite^] 
lattice. Vq, K, e, and / are the fit parameters. 



In the analysis we followed the fitting procedure, described in Ref. p6| . Four different 
fit algorithms have been applied to the data: 



• uncorrelated fits with the errors of potential values obtained on the original 

14 The lattice sums have been computed numerically on 1024 3 , 2048 3 , and 4096 3 lattices and ex- 
trapolated in 1/Ls to their infinite volume limits. 
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4.00,17.32 


35/44 
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K 


e 


Rmin,R m ax 


X Z /N DF 


UN 


.523(1)(6) 


.01451( 4)(14) 


.269( 2)(15) 


4.00,41.57 


70/95 


UB 


.519(1)(3) 


.01466( 5)(15) 


.250(10)(15) 


4.00,41.57 


91/95 


CN 


.519(1)(2) 


.01466( 5) (14) 


.251(12)(16) 


4.00,41.57 


92/95 


CB 


.518(7)(8) 


.01467(25)(31) 


.242(48)(62) 


4.00,41.57 


95/95 


/? = 2.74 , L s xL T = 32 4 




V 


K 


e 


Rmin.Rma* 


X z /Ndf 


UN 


.4816(12)(10) 


.00834(6)( 5) 


.217(3)( 5) 


6.00,27.71 


36/51 


UB 


.4816( 9) (27) 


.00834(4)(12) 


.217(4) (14) 


6.00,27.71 


47/51 


CN 


.4816( 9) (31) 


.00834(4)(13) 


.217(4) (17) 


6.00,27.71 


47/51 


CB 


.4817(14)(55) 


.00834(6)(37) 


.218(6)(89) 


6.00,27.71 


48/51 



Table 5: Three-parameter fits according to Eq. JjA. The first column labels the fit algo- 
rithm. In the last two columns, the "best" fit range and corresponding \ 2 /Ndf values 
are stated. The first errors are statistical only, the second errors include systematic 
uncertainties. 



sample (UN), 

• uncorrelated fits with errors calculated for each bootstrap separately with a 
subbootstrap (UB), 

• correlated fits with the covariance matrix computed on the original sample (CN), 

• correlated fits with covariance matrices calculated on each bootstrap separately 
(CB). 



The fit range has been adapted automatically. For each range, a quality parameter, 

N DF K 

Q = a N A? ' ( 45 ) 

JVDF.max Ail 

of the fit has been computed from the confidence level, a. The largest quality cor- 
responds to the "best" fit range. As a systematic error we have taken the scatter 
between the fit parameters from fits with Q > |Qmax- 

In Tabs. |] and || the results are listed for three- and four- parameter fits, respectively, 
together with their statistical and systematic uncertainties. In addition, the "best" fit 
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13 = 2.5 , Ls x Lt = 16 4 
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CN 
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20/30 


CB 


545(11(31 


.0350(3)(4) 
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245(81(271 


1 73 13 86 


20/30 
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J? R 
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UN 


550(21 ( fil 


.0347(2) ( 5) 


.256(5)(14) 
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2 24 25 


45/60 


UB 


547(21(101 


.0350(2)( 9) 


.245(4)(23) 


34(21(141 


3 00 1 7 32 

O . L/VJ . J_ 1 .OA 


37/46 


CN 


547(21 ( 91 


.0350(2)(12) 


.245(5)(24) 


34(21(151 


3 00 17 32 

K) • VJ VJ ■ -L 1 . ' 7 — 


37/46 


CB 


549(31(141 


.0348(3)(12) 


.251(9)(36) 


33(41(141 


3 00 1 7 1 5 


37/45 






(3 = 2.635 
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X l-LiDF 


UN 


51 9(1 1(21 


.01467(3)(14) 


.254(2) (9) 


29(21(1 1 1 


Q 00 84 


67 /Q5 


UB 


521(11(21 


.01459(2)( 8) 


.263(2)(6) 


23(31 ( 91 


4.24,41.57 


74/93 


CN 


521(11(21 


.01458(2)( 8) 


.263(1)(7) 


15(21( 81 


4 00 36 00 


71/90 


CB 


523(11(31 


.01451(7)(13) 


.267(4)(9) 


16(71(151 


4 00 27 00 


56/76 






P = 2.74 , L s xL 7 


— 32 4 

OA 










K 


e 


/ 




X l 1 ^Db 


UN 


.4823(4) ( 7) 


.00832(3)(5) 


.2220( 9) (21) 


.21(1)(5) 


2.45,26.00 


38/62 


UB 


.4828(3) ( 8) 


.00830(2)(4) 


.2235( 9) (26) 


.19(1)(3) 


2.83,27.71 


51/62 


CN 


.4828(3)( 5) 


.00830(2)(3) 


.2235( 9) (15) 


.19(1)(2) 


2.83,27.71 


51/62 


CB 


.4827(5)(10) 


.00830(3)(6) 


.2234(14)(34) 


■20(2)(5) 


2.83,26.00 


51/61 



Table 6: Four-parameter fits according to Eq. The first column labels the fit algo- 
rithm. In the last two columns, the "best" fit range and corresponding \ 2 /Ndf values 
are stated. The first errors are statistical only, the second errors include systematic 
uncertainties. 

ranges and corresponding x 2 values are included. Correlated and uncorrelated fits show 
little difference in \ 2 values, which gives evidence that correlations between potential 
data at different R- values are small. In the further analysis we use the results of the 
CN four- parameter fits for the string tension K and self energy Vq. 

As expected from perturbation theory the self energy, Vq, decreases slightly with f3 from 
Vq = .547 at f3 = 2.5 down to Vq = .483 at (3 — 2.7 '4. On the two large lattice volumes, 
the Coulomb coefficients, e, are in agreement with the value 7r/12 w .262, as expected 
in the string picture. On the smaller lattice volumes they come out somewhat smaller. 
This might be a result of the different fit ranges: on the large lattices the fit result is 
dominated by the potential values from separations where the string picture is expected 
to be applicable. The parameter /, used to account for the lattice symmetry turns 
out to be of approximately the same size as e in all cases, indicating that violations 
of rotational symmetry can be understood in terms of the lattice one gluon exchange, 
contrary to SU(3) where values / ~ .6e have been found [0. 
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p 


Vk 


(□) 


- B 2loop(P) 


-B~\(3) 


2.3 


.3690(30) 


.39746 (1) 


5.77 


7.20(13) 


2.4 


.2660(20) 


.36352 (1) 


6.04 


7.28( 9) 


2.5 


.1870(10) 


.34802 (1) 


6.31 


7.36( 8) 


2.5115 


.1836(13) 


.34564 (1) 


6.34 


7.38( 8) 


2.635 


.1208( 1) 


.324308(2) 


6.67 


7.51( 6) 


2.74 


.0911( 2) 


.308721(2) 


6.95 


7.63( 5) 


2.85 


.0662( 6) 




7.25 


7.81( 4) 



Table 7: The "(3 -function" B~ 1 {f3) = <9(lna)/<9(ln/3), obtained by use of the perturba- 
tive two loop approximation, and by the interpolation procedure, described in the text. 
The (3 = 2.3 and (3 = 2.4 values are taken from Ref. f$% p^/ , the (3 = 2.5115 value is 
from Ref. j^J , the [3 = 2.85 value from Ref. ffij. 



4.1.3 Determination of the /3-function 



From Eq. 99 it is evident that the action density diverges like 



function of a (see Eq. |39D, where 

B (9) = -7T ( 4? ) 
am a 

denotes the Callan-Symanzik /5-function. For this reason, we shall set out in this 
section to determine the /3-function within the g 2 = 2N//3 region covered by our 
simulations. 

The /5-function can be expanded in terms of the coupling: 

B(g) = -b g 3 - b l9 5 (48) 

with b = lliV/(487r 2 ) and b x = 34A^ 2 /(3(16tt 2 ) 2 ). From Eqs. g7| and || one finds the 
familiar formula 

a = A£V(/3) (l + O^- 1 )) with f(J3) = e-^ iNbo ((3/2Nb ) b ^ 2b o . (49) 



In Tab. 0, results for the square root of the string tension, \fK, and the plaquette 
expectation value are collected for various j3. The results are taken from the present 
simulation and Refs. |27|, ||, ||, |. From Eq. g| we obtain f(f3)VK = ^A^a) by 
using a 2 = K/k. In Fig. [7| the a dependence of A^r(a) = 19.82A L (a) is shown. As 
can be seen from the non- vanishing slope, within the (3 region accessible by present 
computers, higher order contributions to the asymptotic formula Eq. are important. 
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1 I I I I I I I 

0.05 0.1 0.15 0.2 0.25 0.3 

a sqrt(K) 



Figure 7: A^^{o)^/k versus the lattice spacing a^/s, both measured in units of the 
string tension. The estimates for the A values have been obtained from the perturbative 
two loop formula by use of the bare coupling (upper values) and the (3e scheme (lower 
values). Linear fits to the data are indicated. 



A cut-off parameter A E (a) from the effective coupling /3 E = |(1 — (□)) 1 , introduced 
in Ref. [ 30 1 , has also been computed. This cut-off parameter is translated into the MS 



scheme by the relation A 



11.51A£(a). As can be seen, the slope is substantially 



reduced but asymptotic scaling remains violated. 



The difference between the A -1 (a) sets indicates the size of higher order (perturbative 
and non-perturbative) contributions to the /3-function. Within the present (3 range 
we find the points to fall quite well on straight lines (see Fig. so we use a linear 
interpolation of our data in the region 2.4 < (3 < 2.85, according to the parametrization 

A L \a) = A L \0)+r ] a . (50) 

From the fitted slope, r], we obtain the relation 



and arrive at 

c 



15 However, for small a, the leading order correction should be proportional to g 2 oc 1/ In a, instead. 
So, it is not surprising that the linear effective parametrizations do not extrapolate to the same 
continuum limit. 



24 



The resulting values for —B^ 1 , obtained in the two loop approximation and by the 
above fit, are displayed in the fourth and fifth column of Tab. [7], respectively. De- 
pending on how many of the points we include into our fit, we obtain values 52.2 < 
r)l s/k < 59.5. This systematic uncertainty is incorporated into the errors of the B^ 1 
values in the last column. Watch the difference between the "measured" /3-function 
and the corresponding perturbative expression decrease with the lattice spacing, as it 
should be! 



4.2 Colour field distributions 
4.2.1 General features 




Figure 8: The energy density distribution at [3 = 2.5, R = 8 (r ~ .7 fin) in units of 
the string tension. 

We are now in the position to present a survey on the flux distributions and watch the 
formation of flux tubes with increasing distance between the static sources. 

In Figs. H and |9|, we display the situation at /3 = 2.5 and Ra = 8a ~ .7 fm, for the 
energy and action densities, respectively, in units of the string tension. Notice, that the 
mesh is not equidistant in the perpendicular direction because the off-axis separations 
n± oc (1,1) are included. We confirm the earlier observation j|] that magnetic and 
electric field energies are of similar size (within 20%), i.e. definitely dominated by 
higher order contributions in g 2 . This results in a small energy density in the middle 
of the flux tube, which in the case of Fig. |8] is nevertheless well above noise. The 
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Figure 9: The action density distribution at (3 = 2.5, R = 8 (r ~ .7 fm) in units of 
the string tension. 

vertical axis of Fig. || is expanded by a factor —B^ 1 = — d In a/d In (3 ~ 7.36, relative 
to Fig. [| This is suggested from the form of the sum rules (see Eqs. |39] and 41 



The figures show that this is indeed a reasonable choice. The electric flux tube looks 
distinctly broader around the sources: contrary to the values in the physical region, 
the (self-interaction) values at the peaks of the action and energy density distributions 
roughly equal each other. This observation, which is in accord with the sum rule 
prediction, causes the above optical impression. 



Fig. [10] illustrates the action density distribution at equal physical geometry as Fig. [5] 
but with finer lattice spacing (R = 12 on the j3 = 2.635 lattice). The vertical scale 
of Fig. |Tt] is contracted relative to Fig. |S] by the ratio of the corresponding two B 
values. We observe nice scaling of the fields outside the peaks. The latter diverge by 
the expected additional factor a\^ja\ mh ~ ^" 



The elongation of the flux distribution into a tube is traced in Figs. || 11, and 12 for 
the action density. The physical source separations correspond to .7, 1, 1.35, and 1.7 
fm, respectively. Our data representation avoids use of a smoothing procedure, as has 



been applied in previous work [|3T], |32| . We are in the position to judge the significance 
of the actual data from its intrinsic fluctuations. In this way, we are less bound to 
be deceived by the beauty of some graphic interpolating algorithm. Indeed, given the 
quality of our data, we can follow — for the first time in a lattice simulation — the 
flux tube along distances of up to 1.7 fm or 30 lattice sites at j3 = 2.635 (Fig. [P2|) ! 
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Figure 10: The action density distribution at f3 = 2.635, R = 12 (r « .7 /mj m units 
o/ the string tension. Relative to Fig. ^ the vertical axis has been rescaled by the ratio 
of the corresponding B{(3) values. 



4.2.2 Finite a effects 



In this subsection we will start to discuss the systematic errors on our field mea- 
surements. A prominent effect would be expected from the limitation of the lattice 
resolution a to which we will turn first. 



The comparison shown in Fig. 12 between the action density distribution at a quark 
separation of 1.7 fm, obtained at two different lattice spacings, indicates scaling of 
the results outside the self energy region. The same holds qualitatively true for the 
situation at a distance of .7 fm as can be seen from Figs. ^ and Fig. |10[ Thus, we are 
driven to the conclusion that continuum results can be obtained from quark distances 
as small as eight lattice spacings, at least at positions, separated by more than two 
lattice sites from the sources. 



Let us investigate the situation in some more detail. In Fig. [13] we compare longitudinal 
action flux tube profiles obtained at f3 = 2.5 and j3 = 2.635 for r « .7 fm to each other. 
One source is placed at the origin of the coordinate system. The data are appropriately 
scaled in units of the string tension and in addition divided by the expected anomalous 
dimensions from Tab. [7]. The situation is displayed for x± ~ .17 fm. The latter distance 
corresponds to n± = 2 and n± = 3 for the two (3 values, respectively. Though the data 
are compatible to each other within errors, the values obtained at the finer resolution 
tend to be systematically below the corresponding (3 = 2.5 values. The same is found 
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Figure 11: Same as Fig. |^ for R = 12 (r ~ 1 //nj and R = 16 (V ~ 1.35 /mj. 

in a comparison of transverse profiles obtained on the two data sets at the center plane 
between the sources (Fig. |14|). However, due to the errors on the string tension values 
and the /3-function (Notice, that the vertical axis has been scaled by a factor K 2 \), a 
relative overall scale error between the two data sets of about 8% is expected which 
easily explains systematic deviations from our expectation. 

We conclude that at separations R > 8 continuum action and energy density distribu- 
tions can indeed be observed on the lattice. This conclusion is further supported by 
the fact that cylindrical rotational invariance is restored (within errors) as can be gath- 
ered from Fig. 14 where the values obtained at plane diagonal sites (multiples of \/2) 



neatly interpolate between the values, measured along a lattice axis. As we will see in 
Section |5j] violations of this rotational invariance are encountered for R < 6, even at 
the center plane. However, for sufficiently small lattice resolution, these violations can 
be understood in terms of lattice perturbation theory and eventually corrected to ob- 
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Figure 12: The action density distributions for quark separations R = 20 at (3 = 2.5 
and R = 30 at (3 = 2.635, corresponding to r w 1.7 /m. The z axis of the first plot is 
expanded by the ratio of the two B{0) values in respect to the second plot to account 
for the anomalous dimension. 

tain the corresponding continuum expressions. This is beyond the scope of the present 
paper, where we are mainly interested in non-perturbative large distance effects. 



4.2.3 Finite size effects 



Lattice results for the heavy quark potential and colour flux distributions are subject 
to finite size effects (FSE). The impact of FSE onto (smeared) Wilson loops is twofold. 
The ground state potential V(R) itself might depend on the finite volume, due to the 
infra red cut-off. Contrary to the perturbative expectation, previous lattice studies of 
the confined phase of SU(2) and SU(3) gauge theories [JD|, |20[ show that this effect 
already becomes negligible for lattice extents as small as L$ ~ 1 fm. In the present 



29 



C\J 

< 

CO 

E 

CO 



1 

0.9 
0.8 
0.7 
0.6 
0.5 
0.4 
0.3 
0.2 
0.1 


-0.1 



beta=2.5 h>- 
beta=2.635 h- 



-1 -0.8 -0.6 



-0.4 -0.2 0.2 
n_1 sqrt(K) 



0.4 0.6 0.8 



Figure 13: Comparison of longitudinal action density profiles at r = .7 fm between the 
(3 = 2.5 data (R = 8) and the (3 = 2.635 data (R = 12) at n±a = .17 fm in units of the 
string tension. The vertical axis has been multiplied by B(f3). One source is placed at 
position (ni,n±) = 0. The second source is located at n\\[K ps 1.6 (outside the visible 
range). 



simulation we are able to confirm this observation by comparing the 16 4 and 32 4 
potential data at f3 — 2.5. 

In addition one might worry about the impact of mirror sources, due to the toroidal 
structure of the lattice: if one places sources at the positions and R, the corre- 
sponding state is virtually indistinguishable from a state created by so called mirror 
sources. Thus, in the case of the (self adjoint) fundamental representation of SU(2), 
one expects — in addition to a non-vanishing overlap of the creation operator with 
a QQ state with separation D(0) = R — overlaps with states of internal separation 
D(n) = R + mLs with m; being (not necessarily positive) integer numbers. Let 
us consider for the moment the "perfectly" smeared Wilson loop (with no overlap 
whatsoever to excited states). One would thus anticipate 

^(R,T) = 5> m e" y ( D ( m » T . (53) 



In strong coupling these mirror copy effects are exponentially suppressed as the linear 
size grows in all directions. For a large planar Wilson loop the leading order behaviour 
is: 

W(R, T) = e~ KRT (l + e -K{L s L T -2BT) + . . (54) 
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Figure 14: Comparison of transverse action density profiles at the center plane of two 
sources, separated by r = .7 fm between the (3 = 2.5 data (R = 8) and the (3 = 2.635 
data (R = 12) in units of the string tension. The vertical axis has been multiplied by 
B(J3). 

However, in weak coupling, perturbation theory yields W(R, T) = W(L S — R,L T — T) 
in the non zero momentum sector. At least to the lowest order {0(g 2 /(L|Ly))) the 



zero modes do not obey this behaviour [33]. Their influence might become even more 



important to higher orders, especially in the infra red regime of large R. 

The naive geometrical expectation of Eq. |53| is not borne out by the data, neither 
for the potential nor for the action and energy densities. This can be inferred from 
selection rules due to symmetries of the creation operator. As we shall show in the 
following, this indeed happens in case of the Wilson loop, due to a symmetry under 
transformations by center group elements of SU(2) in the fundamental representation: 

z 2 = 

Let us introduce a non-trivial center transformation to all spatial links that point into 
direction i and cross the hyperplane n« = k + 1/2: 

t\ : Ui(n) — > —Ui(n) for all n with rij = k . (55) 

Obviously, the action is invariant under this transformation since each plaquette cross- 
ing the transformation plane contains two such rotated links. 

Our creation operator = Q(0)U(0 — > R)Qt(R) (Eq. |3|) contains the spatial trans- 
porter {7(0 — > R) that is a combination of various paths connecting the two quarks. 
The smearing algorithm (Eq. O) only permits continuous deformations of the straight 
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Figure 15: Differences between the longitudinal action density profiles measured on 16 4 
and 32 4 lattices at [3 = 2.5 for R = 8 and n±_ = 2. 
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Figure 16: Differences between the transverse action density profiles measured on 16 4 
and 32 4 lattices at (5 = 2.5 for R = 8 and n\ = 0. 
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path. Thus, all paths cross the hyper planes rii = 0, . . . , R4 — 1 an odd number of 
times while all other planes are crossed an even number of times. Therefore, rj^JO) is 
an eigenstate of T % k : 



i r t J , < k < Ri 

TfciR "1 4 , elsewhere [bb) 



As the eigenvalues remain invariant under the Hamiltonian evolution they serve as 
conserved quantum numbers. Consequently, in case of the gauge group SU(2), only 
coefficients c m with even rrii are different from zero in Eq. [53]. Therefore, the effective 
"periodicity''^ is 2L$ rather than Ls- For example, the leading order "pollution" for 
the on-axis separation R carries the decay constant V(2L$ — R)- For a linearly rising 
long range potential, these large exponents cause a strong suppression of fake states. 
This explains why such effects remain unseen in the present simulation, even at R as 
large as 3/4Lg. 

The above arguments can be generalized to the local field strength measurement oper- 
ators On( n ) = e i?( n ), 0\R( n )- Or has no overlap to a QQ state, separated by Ls — R. 
The only relevant finite size effects stem from the periodicity 

O r (ji) = O r (L s - m, L s - n 2 , L s - n 3 ) . (57) 

For n taken along the QQ axis, energy and action densities are strongly suppressed 
outside the sources: in case of a dipole field (the leading order perturbative expectation, 
Eq. [52]) the action and energy densities fall off like (|ni| — R/2)~ A for \ni\ > R/2. 
Thus, FSE into the longitudinal direction are negligible. The n± distributions are 
more sensitive to FSE as will be explained in Appendix (see also Fig. |3|). 

A comparison of the potential computed on 16 4 and 32 4 lattices at /3 = 2.5 shows no 
statistically significant bias due to the volume^. The same holds true for the action 
and energy density distributions as a comparison between the two lattice volumes 
shows for the largest QQ separation realized on the smaller lattice, R = 8 (where FSE 
should be strongest). As an example, the two data sets are displayed in Fig. [15] for the 
longitudinal slice n± = 2. Fig. |16| shows the corresponding transverse distributions for 

Tli = 0. 

Since the (3 = 2.635 lattice is comparable in physical size to the 32 4 lattice at (3 = 2.5 
while the (3 = 2.74 lattice has about the same physical extent as the 16 4 lattice, we 
conclude that all our lattices are sufficiently large for the present purpose and that 
FSE are below the statistical accuracy of the present investigation. 
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Figure 17: Comparison between smeared and unsmeared finite T energy density ap- 
proximants in the center (n% = n± = 0) between two sources, separated by R = 4 at 
(3 = 2.74. The horizontal axis is the time T. 

4.2.4 T-stability 



In the two preceding subsections limitations in the lattice geometry have been discussed 
to substantiate the relevance of our lattice results to continuum physics. Here, we 
address the reliability of our ground state results in view of the (necessarily) limited 
temporal extent of the lattice operators. We will explain in some detail how the 
(T — > oo) results shown above have been obtained. 

For the potential from unsmeared Wilson loops, one has to take T ^> R in order to 
obtain asymptotic results as a consequence of Eq. [18|. In the case of field strength 
operators this amounts to T > 2R since excitations are damped by a decay constant 
AV(R)/2 only (instead of AV(R)). For large R values, in which we are interested, it is 
practically impossible to obtain signals at sufficiently large T. However, this situation 



is considerably improved by the smearing procedure, described in Section |2.3| . This 
is evident from Fig. [17| that presents a comparison between smeared and unsmeared 
results for the energy density in the center between two sources, separated by R = 4 



16 Obviously, the coefficients c m can differ from each other, depending on the path combination, 
appearing in the transporter U(0 — > R), unless Ri — mLs- 

17 The ground state overlaps tend to be smaller on the larger lattice, though the same smearing 
procedure has been applied. We suspect that the number of smearing steps needs to be increased 
when working on larger lattices due to a more extended wave function. However, within our statistical 
accuracy, this is not in the colour field distributions. 
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at P = 2.74. Notice the logarithmic scale! 



To all our (smeared) data we have performed four parameter fits according to Eq. ^5] 
as well as two- and three-parameter fits of the form 

(n(S)) w = <n>|o.j2>-|o> 

+ Cl e- nT/R cosh(2vr S/R) (58) 
+ C2 e-^ T -^)/ R , 

where (n)^,^)-^) an d Q are the fit parameters. In case of two parameter fits, c 2 
is constrained to zero. We note that, because of the different temporal positions of 
magnetic and electric insertions (i.e. different values of S at fixed T), the fits have 
been performed separately before combining the expectation values for £ and B to the 
energy and action densities, e and a. 

In all cases, the agreement with our data was remarkable with x 2 /^df values close to 
one. For the two parameter fits we had to exclude the T = 1 data point. The best 
results have been obtained with the three-parameter fits. In case of four parameters, C3 
was found to agree with zero within the (large) statistical uncertainty. Within errors, 
the T — > 00 extrapolated values coincided with the T = 3 value for large R and the 
T = 4 value for small R in all cases. All our results refer to the extrapolated values 
whose errors have been obtained by the bootstrap method 



In Figs. and [19| we exemplify the time dependence of the electric and magnetic 
energy density estimates, £ and B, at (3 — 2.5 for a quark separation R = 6 at the 
position n\ = 0, rij_ = 3. The corresponding two-, three-, and four-parameter fits are 
included, together with the T — > 00 extrapolated values. Due to the early ground state 
dominance, the fits yield fairly stable results. Notice, that due to the fact that the 
distance S of the plaquette insertions from the central time slice alternates with T, the 
parametrizations are discontinuous. For this reason, the fit values are just indicated 
at integer values of T. In case of integrated quantities, needed for computation of 
the width of the flux tube and comparison with the energy and action sum rules, the 
summation was first performed over the electric and magnetic energy densities for fixed 
T, separately, and the T-extrapolation was carried out subsequently, before combining 
the components to the energy and action densities. 



5 Physics analysis 



Having presented and substantiated our numerical results we are now ready to enter 
the physics analysis. 
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Figure 18: Finite T approximants to the electric plaquette expectation value Sq(0,3) 
in presence of two quarks (diamonds), separated by R = 6 at (3 = 2.5 at the position 
Hi — 0, n_i_ = 3 in lattice units. Two-, three- and four-parameter fits are indicated, 
together with the corresponding extrapolated asymptotic values (rightmost points with 
error bars). 



0.0012 



0.001 - 



0.0008 - 



0.0006 



0.0004 



0.0002 




1 



ex 



Figure 19: Same as Fig. but for the magnetic plaquette expectation, £> 6 (0,3). 
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Figure 20: The central transverse energy density profile at R = 2, (3 = 2.5 together 
with fit curves of methods (1) and (4). 
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Table 8: Results of fits to the central transverse profile of the energy density distribution 
at f3 = 2.5. ci, c 2 , and 5 are fit parameters. A is the integrated area. 
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Figure 21: The central transverse energy density profile at R = 4, (3 = 2.5 together 
with a fit to the lattice expression, fi (method (1)). 
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Table 9: Same as Tab. for the energy density profile at (3 = 2.635. 
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Table 10: Same as Tab. H/or the action density profile at (3 = 2.5. 
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Table 11: Same as Tab. $for the action density profile at (3 = 2.635. 
5.1 Transverse shape 



We will focus on the transverse profile of field distributions in the center plane of the 
flux tube. For small separation of the sources, r, perturbation theory is likely to apply, 
and one might thus expect the (energy and action) distributions to follow the shape 



of the dipole field (see Eq. [32] ) 

fd(n±) 



45 2 



(59) 



(47r)(^ + ni) 3 ' 

where the width of the flux tube would increase linearly with R: 5 = R/2. In lowest 
order this would be multiplied by Cpa with a = g 2 /(4fr). 

For small R, the continuum form of Eq. |59] needs be replaced by a lattice sum, fi(n±), 
that can be computed from Eqs. 1531 and 153. Remember, that this is only the leading 
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Figure 22: Same as Fig. 21 at (3 = 2.635 



order perturbative expectation. The data reveals that restoration of rotational invari- 
ance takes place at unexpectedly small separations, especially in the action density. 
To account for these higher order effects, that cancel lattice artefacts in a subtle way, 
we will also allow for a mixture of both, lattice and continuum expressions. 

As the source separation becomes large, compared to the transverse size of the object, 
the string picture comes into play and we might expect (at least for small n±) the flux 
distributions to be proportional to 



1 



ni 



2^ 6XP 



The normalization has been chosen such that 



J2flM ~ f d2n ±fcM = Jd 2 n ± f g (n ± ) = — 



(60) 



(61) 



to allow for a direct comparison of the fitted coefficients. The question arises how 
the lattice data might connect between the two regimes. We will attempt to model 
the transition region by fits to the dipole parametrization, f c (n±), with 5 treated as 
free parameter. This is motivated by the idea that due to antiscreening of the colour 
sources, their effective charge increases when viewed at increasing distance from the 
QQ axis, nx, which is tantamount to a rescaling of R. 



In this heuristic spirit, a wide variety of one-, two- and three-parameter fits (with free 
parameters C\,C2, and 5) have been performed on the data, which are listed here in 
shorthand notation: 
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Figure 23: The central transverse action density profile at R = 8, (3 = 2.5 together 
with fit curves of methods (4) and (6) (dipole and Gauss). 

1- c 2 fi(n ± ), 

2. ci/d(n±) with 5 = R/2, 

3. cif d (n±) + c 2 fi(n±) with 5 = R/2, 

4. ci/ d (n_L;5), 

5. Ci/ d (n_L; 5) + c 2 fi(n_L), 

6. Cif g (n±] 5), and 

7- Ci/ 9 (n_L;5) + c 2 /;0-l)- 

The stable fit results are collected in Tabs. which also contain the integrated 

area, 



25 2 \ R 2 

The formula is only exact in the infinite volume limit. The second term has been 
corrected by numerical computations of the corresponding lattice sums. In case of the 
Gaussian and dipole parametrizations, additional fits, according to the finite volume 
expressions, derived in Appendix (Eq. |118 ), have been performed. Subsequently, the 



results have been corrected for the finite volume in the way, described in the appendix. 
For the Gaussian profile the finite size corrections on the integrated area are negligible 
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Figure 24: Same as Fig. [|3| for R = 12. 

(below .2%). However, in case of a dipole distribution, though 5 is little affected by 
FSE (up to 4%), the impact on the area is substantial (up to 25% at large R !). 

In case of the dipole fits to the energy density, the combination c\ + C2 = Cj?o; c (i?) 
amounts to a kind of effective coupling on the scale R. Notice, that the odd numbered 
ansatze incorporate the lattice expression, fi, while the forms (2), (4), and (6) only 
involve continuum formulae. Ansatze (1) and (2) require one parameter only, while 
(3), (4), (6) are based on two and (5), (7) on three parameters. 

Energy profile We start with the discussion of the energy density data. We will 
concentrate the analysis mainly on the preformation of flux tubes, along the guidelines 
of perturbative prejudice. For our energy density resolution is not yet high enough to 
map out the proper string region, r > .75 fm. 

At (5 = 2.5 and R = 2 the data is very precise and excludes the one-parameter fits 
(1) and (2) as well as the two parameter fit with constrained width (3). The first 
acceptable results are reached with ansatz (4), yielding a widthQ 5 ~ .9 The situation 
is visualized in Fig. |20] where we compare ansatze (1) and (4) against the data^]. 



18 The deviation from the expected value, 5—1, can be attributed to the fact that for small R the 
lattice dipole tends to be more narrow than its continuum counter part (Fig. ^) . 

19 Excluding the point n± = 0, we also find acceptable fits with methods (2) and (3). The same 
is the case at j3 = 2.635 where, due to the link integration procedure, no data point is available at 
this position. It is interesting to see from the large coefficient c\, that the data prefers the continuum 
dipole over the lattice dipole. 
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This situation changes at R = 4 where ansatz (1) leads to good results at both (3 values 
(while (2) fails), as can be seen from Figs. |21"| and [22]. Ansatze (3) and (5) yield results 
of equal quality with c\ <C c 2 which fits very nicely into the lattice perturbative picture. 
The data can also be parametrized by a continuum dipole with width 5 ~ 1.55 (ansatz 
(4)). However, the result of ansatz (5) (ci <C c 2 ) shows that the data prefers the 
(one-parameter) lattice expression to the (two-parameter) continuum expression. At 
R = 6 statistical errors allow for all parametrizations (apart from occasional numerical 
instabilities). 

We conclude that qualitatively the R = 2 data is described by leading order lattice 
perturbation theory while the R = 4 and R = 6 data can be quantitatively understood 
along this line. 

It is gratifying to see that the effective coupling parameter a c (R) increases with R, 
as is expected from asymptotic freedom. For (3 = 2.5 we obtain values ranging from 
.06 < a 2 (2a) < .075 (under exclusion of = 0) over .105 < a 2 {Aa) < .125 up to 
.145 < a 2 c {6a) < .225 while at = 2.635 we find the ranges .055 < a 2 c (2a) < .060, 
.065 < a 2 c (Aa) < .085 and .085 < a 2 c (Qa) < .13, respectively. We note, that 4a 2 .5 ~ 
6a2.635- Thus, these numbers give a consistent picture and should be put in perspective 
to the bare couplings a = 2N/(4^/3) « .127 and a « .121 for (3 = 2.5 and (3 = 2.635, 
respectively. 

Action profile In case of the action density, a pure lattice Coulomb ansatz is expected 
to fail since the action density is largely due to higher order effects. Nonetheless, it 
would be interesting to see whether an admixture of this term within the parametriza- 
tion remains necessary to account for lattice artefacts. 

It turns out that this heuristic approach looks little promising as all fits to the R = 2 
and R = 4 action data yield values x 2 ^ -Wd.fQ Among the fits, the three-parameter 
forms (5) and (7) come closest to being successful. The fits are not good enough, 
however, to decide whether this gives genuine evidence for perturbative lattice effects 
or trivially reflects the higher flexibility of a three-parameter ansatz. 

From R = 6 up to R = 10 (R = 8 at (3 = 2.635) the dipole fits with unconstrained 
width (4) appear to be the best parametrization of the data. Beyond these i?-values, 
we observe the data being equally well described by ansatz (4) and (6) at /3 = 2.5, 
while at /3 — 2.635 from R = 10 onwards the Gaussian parametrization turns out to 
be more robust than the dipole ansatz against statistical fluctuations. From R = 12 
onwards the other fitting methods also started yielding \ 2 w N^p values. Since these 
fits are unphysical in this region we have not included them into the table. 

The quality of the dipole (4) and Gauss (6) fits is exhibited for source separations 
R = 8 and R — 12 at (3 = 2.5 in Figs. ^3] and |24], respectively. 

20 Thc (3 = 2.635, R = 2 data is exceptional since in this case we have omitted the n± = point 
from our fits. 
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In conclusion, leading order lattice perturbation theory is found to describe the energy 
density data well at small R. The fitted amplitude is in accord with asymptotic 
freedom. The lattice dipole term helps finding a parametrization for action density 
data at small R, though it is not a dominant term. Continuum dipole fits to the 
action yield acceptable results from distances of about 0.5 fm onwards. Up to 1 fm 
this continuum parametrization has a width larger than R/2. This effect is at variance 
with the antiscreening picture of colour sources and might well be a lattice artefact 
since a lattice dipole is broader than a continuum dipole in this R region (see Fig. [|). 
For larger R, the combination 25 /R decreases to values substantially below 1. From a 
separation of 1 fm onwards the Gaussian parametrization yields an equally good (and 
occasionally superior) description of the data. 



5.2 String formation 

When the sources are adiabatically pulled apart, the accretion of action density, Aer, 
should — according to the string picture — be strictly localized in the center plane 
between the sources. This holds for R large enough compared to the other inherent 
length scales in the problem, i.e. the transverse width of the tube and the size of the 
Coulomb dominated region. It is not a priori obvious when this .R-asymptotia sets 
in. The accretion phenomenon will be exploited to determine this transition point to 
genuine string formation. 

For this purpose, we differentiate the action density distributions with respect to an 
increase in the source separation. This is done by computing the change, Aor = 
°~r ~ a R-2, under stretches (R — 2) — ► R. 

In Fig. ^ we display the results for (3 = 2.5 and R = 6, 8, and 10, respectively. At 
r = 10a « .85 fm, we find impressive evidence that Aa is in fact zero outside the 
center plane! This does not hold at smaller separations, where Act exhibits a net flow 
of action into the center plane from the next neighbour planes. This latter feature 
is in accord with the dipole picture described in Section |3.1| . This action flow is a 
substantial effect at R — 6 and decreases to the 5% level at R — 8. Within our 
resolution we thus conclude, that the transition point to string formation is located 
atQ R w 9. 

We might view the above analysis as a differential diagnosis of flux tube formation, 
which provides much more sensitivity than the more conventional global tool offered 
by the potential. In fact, the R = 9 transition point into the string regime appears 
to be rather deep in the asymptotia of the linear part of the confining potential: it 
corresponds to the point R^/K = 1.7 in Fig. ||, which in physical units is .75 fm! 



21 Strictly speaking, there is of course no transition point into the asymptotic regime, since the 
transition is smooth! 
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Figure 25: Differences, Aa R , between the action density distributions at R and R—2 for 
R = 6,8, and 10 (j3 = 2.5). The sources have been aligned. The labelling corresponds 
to lattice units. 
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5.3 Sum rules 
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Figure 26: The QQ force, obtained from the potential at [3 = 2.5, compared to the 
integrated center plane action density, scaled by the anomalous dimension B in lattice 
units. In addition to the numerically integrated data (num), the dipole and Gauss fit 
results from Section \5.1\ are displayed. 



Within the string regime, Ra > .75 fm, one can write a differential version of the sum 
rules: 



1 



F(VW^l) = ^(V(R-l)-V(R+l)) 



(63) 



-- E (e R . 1 (0,n 2 ,n 3 ) + e R+1 (0,n 2 ,n 3 )) + O((3- 1 ) (64) 

712,113 



-air / dxx (e.R_i(0, x) + e^ + i(0, x)) 
Jo 



(65) 



For the second equality we have assumed rotational invariance and large x cu t- Likewise 
we can write a differential action sum ruleQ: 



F(VR 2 - 1) w airBifi) / dxx (a R -iifd, x) + a R+1 (0, x)) . (66) 

Jo 



Notice, that the fundamental assumption made for the differential version of the sum 
rules is only justified for i? + l>9at/5 = 2.5. Also, the data has to exhibit 



22 The anomalous dimension of the action density outside the sources and adjacent sites equals 
B(f3), independent of the position, n. 
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approximate rotational invariance and x cu t has to be chosen sufficiently large. We 
start from numerically integrating the data. In varying x cut we try to find a plateau. 
For R + 1 < 10 a clear plateau is established while for the R+ 1 > 10 data, x cut = 10a, 
the maximal distance for which we have performed measurements, had to be chosen. 



Thus, these values are only lower limits on the r.h.s. of Eqs. pa and 66. 



As can be seen from Fig. |26]the action sum rule is consistent with our data for R+l > 6. 
This lends further support to the asymptotic character of our data (in T). Violations 
of rotational invariance appear to be small beyond the two directions on which we 
have performed our measurements. Consistency of the energy density data with the 
sum rule is found, albeit within reduced statistical accuracy. 

The fitted integrated action density, A (rescaled by the factor B), obtained in Sec- 
tion |5.1| (Tab. [Top , is also shown in Fig. [26] for a Gaussian and a dipole transverse 
shape of the flux tube. The Gauss values are substantially smaller than suggested by 
the force. This discrepancy can only be due to a slower n± — > oo fall-off of the data 
than assumed by the Gauss ansatz. Notice, that the string picture is only applicable 
for small transverse fluctuations while a large portion of the integral stems from the 
area of large n±. The dipole values are consistent with the force for large separations 
R. It remains to be clarified whether this is just a lucky coincidence. At least for 
.75 fm< r < 1.5 fm the large n± data also seems to be underestimated by a dipole 
shape. The functional form of the profile can be studied in more detail by varying the 
transverse spatial volume and exploiting the observed FSE (Appendix [B]). 



5.4 Width of the flux tube 



In addition to the (parametrization dependent) results on the width of the flux tube of 
Section |5J] (Tabs. |8]-[TT|), we attempt to compute this important parameter by direct 
numerical integration: 

2 = ^dn ± njq(0, n ± ) 



The results, including their systematic errors from varying n cut , are displayed in Fig. P7| , 
together with the expectations from the above dipole and Gauss fits. We realize that 
this method is not a viable way to determine the R dependence of 5: the relative error, 
Ad, of the numerical integration is intolerably large and the two fit results also differ 
by a factor of about 1.5. This, of course, is related to the large weight with which 
large n± points contribute to Eq. |67|. For R > 10 the data is well described, both by a 
dipole and by a Gaussian parametrization for our n± window within statistical errors, 
and yet the two parametrizations differ substantially at large n±. 

The data on the numerically integrated widths for physical distances below .5 fm (the 
largest separation at which numerical integration of the energy density data could be 
performed) exhibits that the energy density values fall onto the line 8 = R/2 while the 
action density values are significantly larger. This tendency has also been observed in 
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Figure 27: The width of the flux tube, 5, against R at (3 = 2.5 (in lattice units). The 
dashed line corresponds to 5 = R/2. In addition to the numerically integrated results 
from the energy and action density distributions, fit results to the action density from 
Tab. [77| and results from the geometric method (with "f—l) are displayed. 



Ref. and is consistent with our fit results (Tabs. § and |10|) . 

An alternative approach to study the functional dependence, 5(R), is to constrain the 
center plane analysis to the results of the differential action sum rule. In addition, 
we apply a geometric method that will correlate results from different i?-values to the 
extent that we end up with reduced relative errors and all uncertainty cast into a large 
overall scale error. To quote the assumptions: 



1. Accretion of additional energy and action when pulling the sources apart is 
localized in the center plane. 

2. At sufficiently large R, the change of the transverse shape under variations of R 
can be absorbed into two (independent) scale transformations. 



Assumption (1) has been verified from our data in Section [5]2] for distances r > .75 fm 
and f3 > 2.5. Within our statistical errors, assumption (2) is also fulfilled in this 



region, according to the fit results in Tabs. 10 and pT 



In this case, we can define: 

5 2 = 7 4 ■ (6* 
Tin 
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Figure 28: The width of the flux tube, 5, against R in units of the string tension, 
obtained by use of the geometric method (j — 1) for the energy and action density 
distributions at (3 = 2.5 and (3 = 2.635. The vertical line indicates the lower limit of 
applicability of the geometric method. The dashed- dotted curve is the string picture 
expectation (Eq. \$%) for R ^K = 1/3, 5 ^K^y = 0.5. 

A is the area below the curve. It can be fixed by the sum rules with high accuracy. 
For the action density we obtain A from 



A = B~ L ((3)(V(R)-V(R-1)) 



(69) 



Note, that here we have taken an asymmetric derivative of the potential. In case of 
the energy density, A directly equals the force, up to a renormalization constant, h 
denotes the action/energy density in the middle of the tube (n = 0). 7 is a geometry 
factor. Depending on the parametrization it can take the following values^ 



7 = I for a distribution, constant for n± < n max and zero outside of this circle, 
7 = 1 for a Gaussian shape, 
7 = 2 for a dipole shape and 
7 = 3 for an exp(—c\n±\) shape. 



23 Note, that these values only apply to the infinite volume case. At large S/Ls, they tend to be 
smaller. 
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By employing the definition Eq. ^8], a large portion of the error on the width is cast 
into the (overall) uncertainty of the geometry factor 7. 

In addition to data, obtained by use of the other methods, we have included the data 
from this geometric method into Fig. ^7|for the case 7 = 1 (triangles). The differences 
between these points and the Gauss fits (crosses) reflect the fact that the large n± 
data is not well approximated by the Gaussian form. Remember, that this very effect 
has also led to an underestimation of the force in Fig. BR 



In Fig. [28|, we display our geometric results (7 = 1) for the action and energy densities 
at (3 = 2.5 and (3 = 2.635, scaled in units of the string tension. The dashed vertical line 
denotes the distance .75 fm above which the geometrical approximation is justified. 
As can be seen, the data exhibits scaling even below this limit. The width of the 
energy density starts out to be smaller than the width of the action density, as has 
been observed in Ref. |§, but increases faster. It then reaches the same magnitude as 
the action density width, before it disappears under the noise level at about .8 fm. 

Above r = 1 fm the action density data is in agreement with a constant value 



5 a a « 1.17^/7/k » 0.52 fm</y , (70) 

where we expect 7 to take values between one and two. Logarithmic fits to the r > 1 fm 
data according to the string picture expectation, Eq. |37|, yield values 



R VK < 1/3 . (71) 
No lower limit is imposed on the cut-off since the data is also in agreement with a 



constant. A curve with parametrization RoyK = 1/3 and 5q\JK/^ = 0.5 is indicated 
in Fig. Ej| 

The strong parametrization dependence of the r.m.s. width, 5, is reflected in the dif- 
ference by a factor of about 1.5 between the Gaussian and the dipole parametrizations 
(and by the different geometry factors, 7). The half width, p, is much less sensitive 
to the parametrization: both forms are valid interpolations of the data in the small 
n± region. For the two parametrizations, p can be connected to 5 by the following 
relations: 



P 



2<5v / h72^ 1.675 (Gauss) , (72) 



p = 25V2 1 / 3 - 1 « 1.02 5 (dipole) . (73) 



The resulting half widths for energy and action densities in units of the string tension 
are displayed in Fig. |29|. We have attempted a finite volume correction to the dipole 
results by fitting the data to the functional form, described in Appendix [D] (Eq. |118|) , 
and subsequently converting the resulting 5^ values into p via Eq. |73[ This amounts 
to a reduction of S by less than 10%. Differences between the uncorrected dipole data 
and the Gauss results (up to 6%) reflect the systematic uncertainty due to the form 
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R*sqrt(K) 

Figure 29: The half widths of the energy and action flux tubes at [3 = 2.5 and (3 = 2.635, 
p, against R in units of the string tension. The line corresponds to the small R dipole 
expectation. 

of the interpolating curve. We observe nice scaling between both (3 values. We also 
confirm the width of the energy flux tube to be smaller than the width of the action 
flux tube for distances below .5 fm. Both densities increase till r ~ 1.1 fm. The action 
density saturates at the level p ~ .7 fm. 

We conclude that the data beyond 1 fm is in agreement with a constant but does 
not contradict the expected string picture behaviour either. However, the ultra violet 
cut-off, r-Q 1 , of the effective string theory is found to be larger than or 1.3 GeV. 
This has to be compared with a lattice resolution of 2.35 or 3.64 GeV at (3 = 2.5 and 
(3 = 2.635, respectively. Thus, we are not too far away from the limit of applicability 
of the string theory. 



6 Summary and conclusions 

We have demonstrated that Wilson loop plaquette correlations offer a viable access to 
a lattice study of the flux tube problem on the required length scale of one to two fm. 

Prior to a workable application of this tool one must ascertain an essential improve- 
ment of the lattice observation technique: the crucial ingredient of our method is the 
smearing of the parallel transporter within the bilocal QQ creation operator. This 
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secures a controlled ground state preparation of long flux tubes within few lattice time 
slices. Smearing is combined with integration on the time-like links of the Wilson loop 
to cut noise. 



As a result, we can observe flux formation in the action density over lengths well be- 
yond 1.5 fm with spatial resolution .05 fm. We find that — due to a center group 
symmetry of the Wilson loop — finite size effects remain well below the level of accu- 
racy reached in the present simulation, at least as long as L$ is kept larger than 1.3 fm 
and R < \L S - In particular, there are manifestly no effects of /^-periodic distortions 
of the field distribution or potential due to mirror sources^ This implies that we 
can safely accommodate a flux tube of length 1.9 fm on our largest lattice of volume 
(2.7 fm) 4 ! The energy and action densities exhibit the expected scaling behaviour, and 
are consistent with the potential measurements through Michael's sum rules. 

At small distances the flux tube is corrupted by lattice artefacts, which can be un- 
derstood in terms of lattice perturbation theory. This holds in particular for the self 
energy peak around the sources, whose non-scaling behaviour is well in accord with 
perturbative expectations. 

The transverse r.m.s. width of the action flux distribution in the midplane between 
the sources rises with source separation, r, until it reaches a rather constant level for 
separations between 1 and 2 fm. The physical value for this constant remains model 
dependent and ranges between .5 and .75 fm, as we estimated from a set of transverse 
profiles, supplementing our measurements with various plausible assumptions on the 
large nj_ behaviour. For the half width we find a plateau value of p ~ .7 fm. In the 
preasymptotic domain, the action width is observed to rise by a factor six, distinctly 
majoring the width of the energy distribution, before both reach their (common ?) 
plateau values. A logarithmic increase as suggested by string pictures for the "asymp- 
totic" R region is consistent with our data, suggesting a rather large ultra violet cut-off 
on inverse wavelengths in effective string models, Tq 1 > 1.3 GeV. 

In the range r > .75 fm, we observe a remarkable stability of long flux tubes in the 
sense that field accretion exclusively occurs in the center plane of the tube, as the 
sources are further pulled apart. This is another important quantitative support for 
the flux tube picture. The issue of establishing a definite tube profile (like Gaussian 
transverse shape) will remain a rather elusive subject for any numerical approach like 
ours. 



The present research can readily be generalized to the situation of more than two static 
sources, like three quark sources in SU(3) |35[] or the case of "nuclear chemistry", with 
two quark-antiquark pairs in SU{2). The latter has been studied recently by Green 
and coworkers |36| in the context of hadronic potentials, while the former will help to 
answer interesting questions related to the three-body character of colour forces in the 

24 It goes without saying that an approach based on the measurement of Polyakov lines would 
neither be amenable to such an improvement programme towards "early T asymptotics" nor would 
it be safe from L^-periodic effects from mirror sources. 
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proton. Work along this line is in progress. 



The methods described here should also be useful in the quantitative studies of the 



confinement mechanism in the maximal Abelian gauge [13, 14, 15 
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During completion of this work, we have received a preprint of Haymaker et al. 
They work at (3- values up to 2.5 and refrain from applying ground state enhancement 
techniques. Instead, they attempt T-extrapolations on derived flux tube properties. 
This enforces a smaller R range and implies less control on systematic effects. 
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Appendix 



A Weak coupling expansion of field distributions 



In this appendix, we recall the one gluon exchange approximation to the lattice po- 
tential and compute the leading order perturbative contribution to the electric energy 
distribution. 

With the lattice gluon propagator in Feynman-t'Hooft gauge, 

a,/i-^b,v : ° b ^ , k^ = 2sink^/2 , (74) 
a weak coupling expansion of the Wilson loop yields 

W(R, T) = 1 + C F g 2 ]T (G(R, r' - r) - G(0, r' - r)) , (75) 

t,t'=0 

where just terms, extensive in T, have been kept and the leading term one is the 
expectation value of the loop with the interaction switched off. Note, that we have 
neglected the zero momentum contribution in the calculation that is suppressed by a 
factor 1/(L 3 s Lt). The colour factor Cp can be calculated by contracting the colour 
indices of the SU (N) generators T a (a — 1, . . . , N 2 — 1): 

C F = ^Tr{T a T b }6 ab = ^± . (76) 



Fourier transforming Eq. |74| yields 

pikn 

Gin) = E F , (77) 

k^O 

, 27r L s L s 

L s 11 

2tt L t L t 

= J^ mi ' mA = — 2~ ' " ' ~2~ 

for the real space gluon propagator on a finite lattice. With 

ikR 

G L (R) = EG(R,r) = £— * (78) 

and 7(R) = - lim^oo In (W(R, T)) /T one obtains 

V(R) = -C F g 2 (G L {R) - G L (0)) . (79) 
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By construction, the weak coupling expansion of the quantity (Eq. |23|) , 

w^GSSr 1 ) ■ (80) 

involves only interactions between the plaquette and the Wilson loop. All self-interactions 
of the plaquette and Wilson loop are cancelled by the denominator. 



An expansion of the plaquette yields 

(n) = i- Cl9 2 + --- (81) 

with c\ = 2Cfg 2 (G(0) — G(l)) = Cpg 2 on symmetric lattices. If we are interested in 
the leading order behaviour only, the plaquette can be approximated by one. However, 
at realistic values of the coupling higher order corrections are large. At ft — 2.5 we find 
for example (□) = .65198. This observation inspired Parisi to formulate a programme 
of mean field improved lattice operators |37j . The idea is to split every lattice operator 
into a part that corresponds to (discretization dependent) fluctuations on the ultra 
violet lattice scale and a physical infra red part. 



More recently, the deviations of lattice results from perturbative expansions in terms 
of the bare lattice coupling parameter have been explained as being due to large 
contributions from tadpole diagrams P8 |. This circumstance has revived the interest 



in mean field or tadpole improved lattice perturbation theory and operators. The 
hope is to suppress ultra violet contaminations by dividing every link in a given lattice 
operator by its Monte Carlo mean field value uq. This is supposed to procure early 
asymptotic scaling and reliable perturbation theory predictions. 



A popular choice of u is the fourth root of the average plaquette. Following this 
procedure, we should divide the expression on the r.h.s. of Eq. ^ by the average 
plaquette. However, in the end, we are interested in the combination /3(D)w only. 
Since the plaquette in the action Sw—ftU has also to be divided by its mean field value, 
ft is replaced by a mean field coupling ftuF = ft(^)- Performing both replacements, 
the (□) contributions cancel. In this spirit, the definition of action and energy densities 
in Section |2.4| represents already in itself a tadpole improved definition. Keeping in 
mind that in the last step our operator will be multiplied by ft, it is justified to neglect 
the multiplicative (□) factor in Eq. |8^ even in a region where g 2 depending deviations 
from one are not small. 



The two loops being disconnected in colour space, only singlets can be exchanged. 
Thus, we expect an exchange of two gluons as the leading order contribution. Techni- 
cally, this can be seen as follows: for computation of the product of two (real) traces, 
both possible relative orientations of the Wilson loop and the plaquette have to be 
averaged over. Thus, exchanges of single (bare or dressed) gluons cancel. The same 
holds for a triple gluon vertex that can only be attached with two legs to one loop 
and with one leg to the other. Due to the Lorentz structure of the propagator Eq. [7^ , 
magnetic plaquettes cannot interact by a direct exchange of gluons with the timelike 
links of the Wilson loop. 
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The colour factor of the two gluon exchange between the disconnected loops turns out 
to be: 

±Tr{T a T b }±Tr{T c T d }5 ac 5 bd = ^Tr{T a T b }5 ab = ^ (82) 

By squaring the one-gluon exchange contribution, dividing the expression by a factor 
two to avoid an overcounting of gluon exchanges, and performing the T integration, 
we obtain (Again, only terms extensive in T have been kept.): 



(U(n) i4 ) w = yi(n) 
C F g A 



(G L (n- ri + e i )-G L (n-r 1 ) (83) 



AN 

-G^n-rs + eO + G^n-rs)) 2 , 
where the sources are placed at the positions r x = ^ei. and r 2 = — f ei. 

After averaging over the two plaquettes used for construction of the electric field 
operator and multiplying by 2/3/a 4 , we end up with 

n 4 /P 2/_xx _ n 2 r yi(n)+yi(n-ei) 

a \&i W)|o,fl)-|o> - 9 (84) 

while 

a 4 (5^) W>-|o> = C% 4 ) • (85) 

Many of the higher order diagrams contribute to B as well as to S. Thus, we would 
expect a partial cancellation of higher order effects in the energy density 

£ fl (n) - B R {n) 
e M n ) = o • ( 86 ) 



From 

J2vi(n) = 2(G L (0)-G L (R)) , (87) 

we obtain 

a 4 Cfl(n) = g 2 C F (G L (0) - G L {R)) + 0(g 4 ) » V(R) . (88) 

n 

Note, that to order g 2 the action density equals the energy density. However, the 
action density is expected to deviate much more from the leading order perturbative 
expectation since higher order electric and magnetic contributions are added and no 
cancellations of diagrams occur. 

Perturbation theory yields (up to a divergent self-energy part) 

v(r) = -C F g^ (89) 

for the continuum potential. The associated electric field is given (up to a colour fac- 
tor) by (? _1 V (v (x + r/2ei) — i>(x — r/2ei)). In the continuum limit the differences in 
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Eq. [83] will be replaced by derivatives, yielding exactly this expression. After squaring 
and expressing the result in lattice units, one obtains: 

^_ in 1 ^f nj-SgR/2 ni + SgR/2 
e R [n) - g C F (4?r)2 2. ^ |n _ eijR/2|3 |n + eii?/2|3 

which is just the continuum limit of Eq. EM. 




B Action sum rule 



In order to derive the action sum rule we start from the definition (W = W(R,T), 
Dirac indices and spatial position are suppressed): 

In Eq. |24|, the spectral decomposition of the argument of the sum for < S < T/2 
has been carried out. For the plaquette position outside the loop, i.e. S > T/2, we 
obtain: 

(n(S)) w = const, x e ~Ei{L T -T)/2 cogh ( Ei ( Lt / 2 - S)) + ■ ■ ■ . (92) 
The signal is suppressed with the temporal distance of the plaquette insertions from 



the Wilson loop, S — T/2, times the mass gap, E\ = m A +a m 3yK, |39| in the 
exponent. 

After summing over all S, we obtain: 

( D )w,2 = ( n )|0,iJ>-|0> 

We have made use of the fact that the off-diagonal (i.e. S dependent) pollutions, inside 
and outside the loop are incomplete geometric series. The summation gives, apart from 
the constant parts, only multiplicative 1 — e~ nAV ( T+1 > or 1 — e - nE n L T-T+i) f ac t rs. The 
constant 6^(11) is the sum of all off-diagonal coefficients from the expansions Eq. p4] 
and Eq. [J2|of (□)w, weighted by corresponding coefficients (1— e _A ^)' 1 or (1— e~ Ei )~ l . 

Together with 

(W)=fvUWe f3u , U= £ U, v {n) , (94) 

n,fj,>u 

and the decomposition of the Wilson loop (Eq. |j) we obtain from Eq. [93|: 
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i|lnW (95) 

1 <9 In |do| 2 _ dV _ dAV K| 2 ^_ AyT 
T dp ~d(3 ~ d(3 Jd 



-e 



Note, that the equalities are exact! Thus, in the expansion of the derivative of the 
Wilson loop the same terms appear as in Eq. |93|. 

A comparison of the 1/T coefficients between the above equation and b of Eq. p3| yields 

b = E W») = • 06) 

From the estimate for ground state overlaps of unsmeared operators Eq. [Tj| we obtain 

BV 

B = -R-^ « /3-V Q it: (97) 

for large R and weak coupling. The monotonous increase of the ground state overlap 
at fixed R with (3 is confirmed in the present simulation. Therefore, B is positive. 
For smeared operators, Vq is replaced by some constant / that is small compared to 
all R~ l , such that the exponential can be expanded and the ground state overlaps 
\do(R)\ 2 exhibit the linear behaviour of Fig. Q / is expected to be proportional to g 2 
to the lowest order such that B pa (3~ l fR. Under the assumption that d\ dominates 
other excited state overlaps, we obtain \d\(R)\ 2 pa fR. Eq. 91 can also serve as a 
definition for colour field measurement operators. However, we have preferred to use 
(□)w instead, due to the better asymptotic behaviour: excited states are suppressed 
by factors proportional to \^Re~ AVT ^ 2 ~ y/Re~^ T instead of R/T only. 

From Eqs. [Zg, |2T], |23], g7| and one obtains: 

E (U^))w,2 ^ £ E \ (S(n) - B(n)) = ^ E ^fl(n) • (98) 

By carefully comparing the coefficients of the expansions one finds many (exact) "sum 
rules" . In the following we list three such examples. 

?"VW - ' (101) 

cr^ denotes the action density distribution of the first excited QQ state without angular 
momentum, a A + is the action density distribution in presence of the lightest glueball 
state. 
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The ground state potential consists of a constant physical part v(R) and a self energy 
contribution Vq which diverges in the continuum limit: 

V(R) = av(R) +V (102) 

By using this decomposition, we obtain from Eq. ^9| the action sum rule 

v-^ o , . d In a . „. 1 dVn 

? aV ^ = -^w-^ ■ (103) 



C Energy sum rule 

The derivation of the energy sum rule, though the more intuitive one, turns out to 
be more complicated. We start from the decomposition of the Wilson action into a 
spatial and a temporal part 

S w = -(3 t U t -(3 s U s ■ (104) 

In the following, the spatial and temporal lattice spacings will be called a s and a t , 
respectively. The asymmetry is defined by £ = a s /a t . Following the steps of the 
previous section, one finds: 

|:£^(n) ^ £E(^(n))w, 2 = |^ln(W) , (105) 

Z P n n ,i 1 °Pt 

|-£^(n) ^ -E W n )>W,2 = -~MW> . (106) 

/ P n n ,i>j 1 °Ps 



A weak coupling expansion [4(J relates the anisotropic lattice couplings to the isotropic 
coupling /3(a s ): 

ftf = (3 + 2Nc s (0 + O(f3- 1 ) , (107) 
M = P + 2Nc t (0 + O(p- 1 ) . (108) 
The coefficients fulfil the relations: 

lliV 

c.(l) = ct{\) = , c ' s (l) + c ;(l) = 6 = — . (109) 



The derivatives of the coefficients have been calculated by Karsch in Ref. | i0| . For 
SU{2) the result is 

c' s = c' s (l) = .11403... , c' t =4(1) = --06759... . (110) 

After expressing the derivatives in respect to the asymmetric couplings by derivatives 
in respect to (3 and £ and taking £ = 1, we end up with: 



60 



Subtracting both expression yields the (exact) action sum rule Eq. ||| 

Adding Eqs. |111| and |112| and dividing by a factor two yields the energy sum rule 



£a 3 e R (n) = ± ( V ( 1 - £(<4 - c' s ) ) + ^^>o + • • •) (113) 




+ Knr^ + 7^h^o + --- (H4) 



The energy sum rule is not exact due to the perturbative origin of the relation be- 
tween the symmetric and asymmetric couplings (Eqs. |107| , |108|) . Of course, it would 



be preferable to measure the corresponding derivatives of In a directly on the lattice 
instead. 



Note that the coefficient of the last term in Eq. |114| is identical to the action sum 



Eq. |103| . The factor d In a/d In (3 appearing in front of v within this term carries an 



anomalous dimension (as the action does), cancelling a factor. Thus, an additional 
factor —v/4 seems to survive the continuum limit a — > 0. Its origin is an incomplete 
resummation of the series: the non-perturbatively determined coefficients contribute 
to all orders in The order (3~ l term yielding the above — v/4 contribution has to 
be cancelled by other anomalous terms appearing in higher orders of the expansion. 
However, if consistently cutting the expansion at order (3~ 2 by expanding the potential 
perturbatively, the coefficient one is reproduced in accord to the expected continuum 
limit. 



D Finite Size Corrections 



In this appendix, we elaborate details on the computation of finite size corrections on 
the action/energy density distributions within the center plane. These FSE are mainly 
due to the periodicity of the measurement operator, 

O R (0,n 2 ,n 3 ) = O R (0,L s - n 2 ,L s - n 3 ) , (115) 

at a given QQ separation, R = Re\. This mirror source effect should not be confused 
with contributions from the replacement R — > R±nLs which are negligible (as shown 
in Section [4.2. 3|) . We also neglect effects from mirror sources along the QQ axis which 
are almost completely screened from the center plane since the colour field densities 
fall off at least as fast as —R/2)~ A into the longitudinal direction. The effect from 
mirror copies placed along the transverse directions can be substantial (depending on 
the ratio 5/Ls). 
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We perform our calculations for two models, namely a dipole transverse shape, 



c <r 

f d (x ± ; oo) = - 3 , (116) 

7T (d 2 + Xj_) 



and a Gaussian shape 



x 2 : 



/^oo) = ^exp • (H7) 



fd(x±] oo) and f g (x±] oo) are the corresponding (infinite volume) center plane en- 
ergy/action density distributions. 

As argued above it is justified to neglect interactions between different pairs of mirror 
sources. We also assume that the chromo electric and magnetic fields on the finite 
volume can be obtained by superimposing the (infinite volume) fields of all (pairs of) 
mirror sources. Note, that we have to add the fields rather than the action and energy 
densities themselves. From the geometry it is clear that the perpendicular electric 
and longitudinal magnetic field components vanish in the center plane. Under the 
assumption that the (perpendicular) magnetic field component is proportional to the 
(longitudinal) electric component, we find: 

f(x 2 , x 3 ; L s ) = ij2 g (x 2 + n 2 L s , x 3 + n 3 L s ) I (118) 

with g(x 2 ,x 3 ) = g(x±) = \J f(x±; oo). The integrated area can be computed in the 
following way: 



A(L S ) = / dx 2 dx 3 f(x 2 ,x 3 ; L s ) (119) 
Jo 

r(m 2 +l)L s Am 3 +1)L S 

= 2^2^ dx 2 dx 3 g(x ± )g(x 2 + n 2 L s ,x 3 + n 3 L s ) 

n 2 n 3 m,2m 3 J ™- 2 Ls Jm 3 L s 

Jd 2 x ± g(x 2 -n 2 L s /2,x 3 -n 3 L s /2) 

xg{x 2 + n 2 L s /2,x 3 + n 3 L s /2) . (120) 



In case of a dipole field with (infinite volume) r.m.s. width, S, we have 

9d(x±)cx (5 2 + x 2 ± y 3/2 (121) 

and obtain for the area, 

r2n rOO . . . . —3/2 

A(L S ) oc d( P rfrr(d 2 (n 2 ,n3)+r 2 (l-L|(n 2 cos0 + n3sin0) 2 )) (122) 

n 2 n 3 J° "'O 

with 

d\n 2 ,n 3 )=5 2 + ^f(n 2 2 + nl) . (123) 
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After performing the radial integration, we arrive at: 

A{Ls) K U 3 h 4* - LUn 2 cos + n 3 sin & (124) 

= - y — 1 (125) 

2 ^ d 3 ^ 2 - L|(n| + n|)/4 

/ r2 \-3/2 

= A(oc)£ U + t|(^ + ^) (126) 



with A(oo) = c/(25 2 ). 
For the Gauss fits we have 



9g(x±) oc exp (^-^j ■ ( 127 ) 



Thus, we end up with 



A(L S ) = E/^exp(-^exp^(n 2 + n^ (128) 



A(oo)]Texp (-^("1 + "!)) 

"2^3 \ / 



In conclusion, the results for both transverse profiles read 



(129) 



ML.) = M°°) T. ^ L >'+"* )/4) and (130) 

n 2 n 3 9d{V) 

AIT \ A(„\ \^ ( L K"l + "2)/2) nql x 

A(L S ) = A(oo) 2^ — — , (131) 

n 2 n 3 9g\y) 



respectively, with g(x±) = y/(xj_;oo) and A(oo) = c/(2<5 2 ). For the typical dipole 
r.m.s. width 5/L$ = 6/32 we find an increase of the area by 43% due to the finite 
volume while the corresponding Gaussian result (5/Lg ~ 4/32) is only affected by 
5 x 10~ 7 . Notice, that the infinite volume S can be obtained by fits of the form Eq. |118| 
from finite volume data. 
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